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Abstract 



We develop Wigner's approach to a dynamical transition state theory in phase 
space in both the classical and quantum mechanical settings. The key to our de- 
velopment is the construction of a normal form for describing the dynamics in the 
neighborhood of a specific type of saddle point that governs the evolution from re- 
actants to products in high dimensional systems. In the classical case this is the 
standard Poincare-Birkhoff normal form. In the quantum case we develop a normal 
form based on the Weyl calculus and an explicit algorithm for computing this quan- 
tum normal form. The classical normal form allows us to discover and compute the 
phase space structures that govern classical reaction dynamics. From this knowledge 
we are able to provide a direct construction of an energy dependent dividing sur- 
face in phase space having the properties that trajectories do not locally "re-cross" 
the surface and the directional flux across the surface is minimal. Using this, we 
are able to give a formula for the directional flux through the dividing surface that 
goes beyond the harmonic approximation. We relate this construction to the flux- 
flux autocorrelation function which is a standard ingredient in the expression for 
the reaction rate in the chemistry community. We also give a classical mechanical 
interpretation of the activated complex as a normally hyperbolic invariant manifold 
(NHIM), and further describe the structure of the NHIM. The quantum normal form 
provides us with an efficient algorithm to compute quantum reaction rates and we 
relate this algorithm to the quantum version of the flux-flux autocorrelation function 
formalism. The significance of the classical phase space structures for the quantum 
mechanics of reactions is elucidated by studying the phase space distribution of scat- 
tering states. The quantum normal form also provides an efficient way of computing 
Gamov-Siegert resonances. We relate these resonances to the lifetimes of the quan- 
tum activated complex. We consider several one, two, and three degree-of-freedom 
systems and show explicitly how calculations of the above quantities can be carried 
out. Our theoretical framework is valid for Hamiltonian systems with an arbitrary 
number of degrees of freedom and we demonstrate that in several situations it gives 
rise to algorithms that are computationally more efficient than existing methods. 
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1 Introduction 

The subject of this paper is transition state theory - classical and quantum. Transi- 
tion state theory (TST) (sometimes also referred to as "activated complex theory" 
or the "theory of absolute reaction rates") is widely regarded as the most important 
theoretical and computational approach to analyzing chemical reactions, both from 



a qualitative and quantitative point of view. The central ideas of TST are so funda- 
mental that in recent years TST has been recognized as a very natural and fruitful 
approach in areas far beyond its origin of conception in chemistry. For example, it 
has been used in atomic physics [JFUOOj . studies of the rearrangements of clusters 
\KB99\ IKBn2j . sohd state and semi-conductor physics |.TTDF84I IEck95j . diffusion 
dynamics in materials |VMG02j . cosmology |dQdAST02j . and celestial mechanics 
[.TRL+n2|IWBWn5b| . 

The literature on TST is vast, which befits the importance, utility, breadth, scope, 
and success of the theory. Searching ISI Web of Knowledge on the phrase "transition 
state theory" yields more than 17,600 hits. Searching Google with the same phrase 
gives more than 41,000,000 hits. There have been numerous reviews of TST, and 
the relatively recent review of |TGK96j is an excellent source for earlier reviews, 
historical accounts, books, pedagogical articles, and handbook chapters dealing with 
TST. Moreover, |TGK96j is notable from the point of view of that in little more that 
10 years it has attracted more than 458 citations (and it also contains 844 references)! 

Certainly the existence of this vast literature begs the question "why does there 
need to be yet another paper on the theoretical foundations of TST, what new 
could it possibly add?" The one word answer to this questions is, "dynamics" . Ad- 
vances in experimental techniques over the past twenty years, such as, e.g., femtosec- 
ond laser spectroscopy, transition state spectroscopy, and single molecule techniques 
( |Neu921 IPZ951 IZewOOj ) now provide us with " real time" dynamical information on 
the progress of a chemical reaction from "reactants" to "products" . At the same time, 
these new experimental techniques, as well as advances in computational capabilities, 
have resulted in a growing realization among chemists of the ubiquity of non-ergodic 
behaviour in complex molecular systems, see, e.g. |SY04l IBHC05t IBHCOGt ICar05] . 
All of these results point to a need to develop a framework for studying and under- 
standing dynamics in high dimensional dynamical systems and recently developed 
tools in computational and applied dynamical systems theory are giving new insights 
and results in the study of the dynamics of molecular systems with three or more de- 
grees of freedom. In particular, we will show how these recent advances in analytical 
and computational techniques can enable us to realize Wigner's dynamical picture 
of transition state theory in phase space for systems with three or more degrees- 
of-freedom. However, to set this in context we first need to describe a bit of the 
historical background and setting of TST. 

Transition state theory was created in the 1930's, with most of the credit being 
given to Eyring, Polanyi, and Wigner, who are referred to as the "founding trinity 
of TST" in Miller's important review on chemical reaction rates ([Mil98b]). Never- 
theless, important contributions were also made by Evans, Farkas, Szilard, Horiuti, 
Pelzer, and Marcelin, and these are described in the discussions of the historical 
development of the subject given in |LK83t [PTOSaj . 

The approach to TST taken by Eyring |Eyr35| emphasized thermodynamics (see 
the perspective article of [PetOO] ) . The approach of Wigner |Wig38| on the other 
hand is based on classical mechanics (see the perspective article of [GarOOj ). It is 
the dynamical approach of Wigner that is the focus of this paper. Despite the fact 
that the original framework of TST is classical mechanics, it is natural to consider 
quantum mechanical versions of this approach to reaction dynamics. We will first 
describe the classical mechanical setting, and then consider the quantum mechanical 
version, and we will emphasize how much of the structure and philosophy of the 
classical approach influences the quantum approach. 



1.1 Transition State Theory: Classical Dynamics 



To begin with, we first examine the assumptions of classical TST, as set out by 
Wigner. Wigner begins by stating that he considers chemical reactions in a set- 
ting where the equilibrium Maxwell-Boltzmann velocity and energy distributions are 
maintained (see |Mah74j for a detailed discussion of this point) and for which the 
potential energy surface is known ( [GarOOj ). He then gives the following assumptions 
from which he derives TST: 

1. the motion of the nuclei occurs on the Born-Oppenheimer potential energy 
surface ("electronic adiabaticity" of the reaction) 

2. classical mechanics adequately describes the motion of the nuclei 

3. there exists a hypersurface in phase space dividing the energy surface into a 
region of reactants and a region of products having the property that all tra- 
jectories that pass from reactants to products must cross this dividing surface 
precisely once. 

It is important to note that Wigner clearly developed his ideas in phase space, the 
arena for dynamics. It is important to keep this in mind since a great deal of later 
developments occur in configuration space, in which certain dynamical properties are 
obscured. 

From the modelling point of view, the first two assumptions are of a very different 
nature than the third. The first two are central to developing the model, or dynamical 
system (i.e. determining the potential energy surface and Hamiltonian function). As 
a result, once a dynamical system describing the reaction has been developed the 
third "assumption" cannot really have the status of an assumption. Rather, such 
a hypersurface satisfying these properties must be shown to exist for the dynamical 
system. Of course, in practice this is exactly how the theory is utilized. One starts 
with a dynamical system describing the reaction, and then one attempts to construct 
a "dividing surface" having the required characteristics. It is precisely this third 
"assumption" that is at the heart of this paper and from which, as we shall see, 
many dynamical consequences flow. 

We will be concerned with dynamics on a fixed energy surface. In this paper "en- 
ergy" means the total energy of the system, e.g the sum of the kinetic and potential 
energies. More mathematically, the energy surface is the level set of the Hamiltonian 
function^ This is important to keep in mind because in not an insignificant portion 
of the relevant literature the meaning of the phrase "energy surface" is actually the 
"potential energy surface" , and a great deal of effort is expended in attempting to 
infer dynamical phenomena from the "topography" of the potential energy surface. 
Certainly for one degree-of-freedom (DoF) Hamiltonian systems (i.e. one configura- 
tion space coordinate and one associated momentum) one can understand all possible 
dynamics from the shape of the potential energy surface. This is definitely not true 
for more than one DoF (or else dynamical phenomena such as "chaos" would have 
been discovered many years earlier). However, two DoFHamiltonian systems, where 
the Hamiltonian is the sum of the kinetic and potential energies, does allow for cer- 
tain constructions based solely on the potential energy surface that imply important 
dynamical phenomena. We will survey these later in this introduction. We empha- 
size that similar constructions using the potential energy surface for systems having 
more than 2 DoFsimply do not work in the same way as they do for 2 DoF. 

^Hamiltonian functions can be more general than the sum of the kinetic and potential energy terms. 
They could contain magnetic terms or Coriolis terms, for example. Nevertheless, we will still refer to the 
level set of the Hamiltonian function as the "energy surface" . 



Now to realize assumption 3, on a fixed energy surface, we need to choose a 
dividing surface that will "separate" the energy surface into two parts ( "two parts" is 
a bit too simplistic, but we will come back to that later) - one part corresponding the 
reactants and the other to products. The dividing surface would have the additional 
(dynamical) property that trajectories evolving from reactants to products cross it 
only once. Again, these reactant and product regions are typically defined via the 
potential energy surface. They are often interpreted as "potential wells" (i.e., local 
minima of the potential energy function) that are "separated" by a "saddle point" 
and a surface (in configuration space) passing through a neighborhood of the saddle 
point serves as the dividing surface ( |Pec81| ). We will show that for systems with 
three or more DoFsuch a configuration space approach, in several different ways, 
does not allow one to realize Wigner's original construction of TST. In fact, this is a 
central message of this paper. It can be misleading, and even wrong, to attempt to 
infer dynamical phenomena from the topography of the potential energy surface. 

In the series of papers 

IWBWOSc] the fundamental framework for phase space TST is developed. The start- 
ing point is classical mechanics and a Hamiltonian function describing the system 
(the same as |Wig38| ). The Hamiltonian can be expressed in any convenient set of 
coordinates, have any number, d, degrees of freedom (DoF), and does not have to be 
of the form "kinetic plus potential energy", e.g., it can include rotational or magnetic 
terms. 

With the Hamiltonian function in hand, the next step is to locate particular 
saddle-like equilibrium points of the associated Hamilton's equations that are of 
a certain type. Namely, the matrix associated with the linearization of Hamil- 
ton's equations about the equilibrium has a pair of real eigenvalues of opposite 
signs (±A) and 2d — 2 purely imaginary eigenvalues occurring in complex conju- 
gate pairs (iiwfc, k = 2, . . . , d). Such equilibria are called saddle-centre-- • • -centres, 
and structures associated with these equilibria provide the fundamental mechanism 
for "transformation" in a large, and diverse, number of applications (some listed at 
the beginning of this introduction), whose dynamical consequences have remained a 
mystery. 

Of course, locating saddles is in the spirit of classical transition state theory, 
but there is an important difference here. We are concerned with the dynamical 
consequences of certain types of saddles of Hamilton's equations in phase space. The 
usual approach is to consider saddles of the potential energy surface (the setting of 
the "landscape paradigm" [Wal04] ) . However, if the Hamiltonian has the form of the 
sum of the kinetic energy and the potential energy, then there is a correspondence 
between the rank one saddles of the potential energy surface and the saddle-centre- 

• • • -centre type equilibria of Hamilton's equations. But here we emphasize the phase 
space setting and the influence of this saddle in the dynamical arena of phase space. 
We reiterate that a central point of ours is that it is difficult, and often misleading, 
to try to infer dynamics from properties of configuration space. 

Next we seek to understand the phase space geometry near this saddle-centre- 

• • • -centre (henceforth referred to as a "saddle" ) equilibrium point of Hamilton's 
equations. An understanding of the geometry will give rise to a set of coordinates 
that will enable us to explicitly compute the phase space structures that govern 
transport and to quantify their influence on trajectories. This set of coordinates is 
realized in an algorithmic manner through the use of the Poincare-Birkhoff normal 
form procedure. These normal form coordinates are central to our theory and the 
resulting analytical and computational techniques. In particular, they enable us to 
show that "near" the saddle the energy surface has what we call the "bottleneck 



property" which facihtates the construction of an energy dependent dividing surface. 
This dividing surface has the "no-recrossing" property and the flux across the dividing 
surface is "minimal" (in a sense that we will make precise). Moreover, the coordinates 
also naturally give rise to a "dynamical reaction path". We want to describe these 
notions in a bit more detail and place them in the context of the chemistry literature. 

Further, we note that historically it has been well-recognized that the computation 
of quantities associated with chemical reactions is greatly facilitated by adopting a 
"good" set of coordinates ( [JRHTI IEM741 IMil76[ lMil77j ). In particular, if the 
Hamiltonian is separable, i.e. there is a set of configuration space coordinates in terms 
of which the equations of motion decouple, then the choice of a dividing surface with 
the no recrossing property is trivial ([ GarOOj ). This situation is extremely special and 
therefore almost irrelevant for chemical reactions. However, the normal form method 
shows that such a decoupling can always be obtained in the neighbourhood of the 
dividing surface through the symplectic ( "canonical" ) transformation of the full phase 
space coordinates (i.e. a symplectic transformation mixing the configuration space 
coordinates and the conjugate momenta). The normal form thus is a constructive 
way of obtaining "good" coordinates in phase space. 

The Bottleneck Property of the Energy Surface and the Energy De- 
pendent Dividing Surface: The geometry or "shape" of a fixed energy surface 
has received little attention, as opposed to consideration of the geometry or "shape" 
of potential energy surfaces. This is unfortunate since an understanding of the geom- 
etry of the energy surface is essential for constraining and interpreting the possible 
global dynamics. Nevertheless, the lack of attention to this issue is understandable 
since such considerations give rise to extremely difficult mathematical problems. As 
an example, the importance of an understanding of the topology of the energy surface 
for an understanding of the dynamics of the three body problem was emphasized by 
Poincare ( [Poi93al IPoi93bt IPoi93c] ) , and work on this problem has involved some of 
the giants of mathematics of the 20*'' century and has resulted in the creation of many 
new areas of mathematical research. Very recent results on the three body problem, 
as well as a discussion of the history of the subject, can be found in |MMW98] . and a 
discussion of the developments of an appropriate computational framework for study- 
ing such questions for general Hamiltonian systems can be found in |KMM04] . We 
would expect that similar studies of the structure of the energy surfaces for standard 
Hamiltonian's arising in studies of reaction dynamics will be similarly fruitful and 
lead to new global dynamical insights. 

However, there are "local" results that describe the geometry of the energy surface 
that are very relevant to studies of reaction dynamics and TST. In particular, for a 
range of energies above that of the saddle, the {2d — l)-dimensional energy surface 
has locally the structure of the product of a {2d — 2)-dimensional sphere with the 
real line, S'^'^~^ x M. We say that in this region of the phase space the energy 
surface has the "bottleneck property" because it is (locally) separated into two pieces: 
^2d-2 ^ ^+ ^^^^ ^2^-2 ^ g^j^j 52^-2 ^ |Q| -g ^Yie dividing surface that separates 

these two pieces of the energy surface, and we identify the two pieces separated by this 
dividing surface as "reactants" and "products" . It will turn out that M corresponds 
to a natural (energy dependent) "reaction coordinate" and 5^"^"^ will correspond 
to (energy dependent) unstable bath modes, or vibrations "normal" to the reaction 
coordinate. 

It should be clear that the geometry of the energy surface, as the energy varies, 
is an important feature of reaction dynamics. In particular, the geometry changes 
with energy and the "bottleneck", S'^'^^^ x M, may deform into a more complicated 



shape as the energy is further increased above that of the saddle. This can lead 
to the "breakdown" of the validity of transition state theory in the sense that we 
are not able to construct a dividing surface separating reactants from products that 
is not recrossed. We note that an "energy limit" for TST has been discussed in 
[(tL77[[SK78] . Looking at it another way, the energy surface deforms in such a way 
that the distinction between reactants and products becomes unclear. This is one 
way in which TST can "break down" . We will mention one other way after we have 
introduced the notion of a normally hyperbolic invariant manifold. 



The "No-Recrossing" Property and Minimal Flux: The dividing surface 
described above can be realized through the normal form computations and trans- 
formations ( [WWJUOll IUJP"'"Ol] ). The high dimensional spherical geometry, S"^"^"^, 
is significant in several ways. A sphere is separated into two parts along its equa- 
tor, which in this high dimensional case is given by S"^"^'^, the {2d — 3)-dimensional 
spherell The Hamiltonian vector field is transverse to each hemisphere, but in an 
opposite sense for each hemisphere. This indicates the evolution from reactants 
to products through one hemisphere, and the evolution from products to reactants 
through the other hemisphere. Transversality of the Hamiltonian vector field to a 
hemisphere is the mathematical property one needs to show that there are "no local 
recrossing of trajectories" , as is shown in UJP^Ol] and in this paper. The Hamilto- 



nian vector field is tangent to the equator of the spherell Mathematically, this is the 
condition for the equator, S'^'^~^, to be an invariant manifold. More precisely, it is 
saddle like in stability and an example of a normally hyperbolic invariant manifold, 
or NHIM ( |Wig94[IWWJU01l[UJP+01p . The NHIM has the physical interpretation 
as the "activated complex" - an unstable super molecule poised between reactants 
and products. 

Except for the equator, S'^'^~^, (which is a normally hyperbolic invariant mani- 
fold), the dividing surface thus is locally a "surface of no return" in the sense that 
all trajectories that start on the dividing surface exit a neighborhood of the dividing 
surface UJP"'"Ol] . Most importantly for reaction dynamics, the energy surface has 



the "bottleneck property". That is, our dividing surface locally divides the energy 
surface into two, disjoint components, which correspond to reactants and products. 
Therefore the only way a trajectory can pass from one of these components of the 
energy surface to the other is to pass through the dividing surface. The issue of "re- 
crossing" is an important part of the choice of the dividing surface. Truhlar [Tru98j 
distinguishes two type of recrossing: local and global recrossing. Local recrossing 
cannot occur with our choice of dividing surface. However, global recrossing is a very 
different matter. If the energy surface is compact (i.e. closed and bounded, for our 
purposes) then the Poincare recurrence theorem ( |Arn78j ) implies that global recross- 
ing must occur for almost all trajectories crossing the dividing surface. Moreover, 
the existence of homoclinic orbits and heteroclinic cycles may also be an intrinsic 
feature of the dynamics QWBW04al IWBW04bl IWBWOSal IWBWOSbl IWBWOScj I. 
Their existence also implies that global recrossing cannot be avoided regardless of the 
choice of transition state; in other words, global recrossing is a fundamental property 
of the dynamics and its presence does not therefore indicate the limitations of any 
particular method for constructing a dividing surface. 



^Think of the familiar, and easily visualizable, case of the 2-dimensional sphere, S^. It is separated 
into two hemispheres by its equator, a sphere of one less dimension, S^. 

■^If the Hamiltonian vector field is transverse to one hemisphere, transverse to the other hemisphe 
in the opposite directional sense, and it varies smoothly in phase space, then we can view the equator 
where the Hamiltonian vector field "changes direction" . 



Wigner Wig38| pointed out that the effect of trajectories recrossing the divid- 
ing surface would result in "too high values of the reaction rate" . This observation 
naturally leads to the notion of variational transition state theory, where the idea is 
to vary the choice of choice of dividing surface in such a way that the flux across 
the dividing surface attains a minimum value (see |Kec67j and the review paper of 
|TG84j ). The latter review paper contains 206 references and has more than 390 
citations, which is indicative of the fact that variational transition state theory is a 
huge subject in its own right. Much of the work that falls under the heading of "vari- 
ational transition state theory" involves dividing surfaces in configuration space (see 
[JJOli IBJOSj for a systematic development of this approach) . These beautiful results 
obtained by our predecessors can with modern day mathematical tools be formulated 
differently, leading to more general results. The beginnings of a general framework 
for such an approach was first given by |Mac91j . and this was used in |WW04j to 
show that the dividing surfaces computed by the normal form approach described 
in [WWJUOH IUJP"'"Ol"] have "minimal flux" . It is worth re-emphasizing again, that 
our dividing surface construction and our flux calculations are carried out in phase 
space, not configuration space. The work in |WW04j implies that one cannot find a 
surface in configuration space for systems with more than 2 DoFthat is free of local 
recrossing, and therefore has minimal flux (unless the system is given in coordinates 
in which the Hamiltonian is separable or has some very special symmetries). 

It is worth pointing out here that flux across a dividing surface is a "local prop- 
erty" with respect to the given surface, i.e. it does not require integration of tra- 
jectories for its computation. If one makes a "bad" choice of dividing surface that 
is not free of local recrossing then one must compute trajectories to correct for the 
local recrossing effect (in the chemistry literature these are referred to as "dynami- 
cal corrections" to the rate, see |MM971 IPriOSj for specific examples of the effect of 
recrossing and how it is treated). This is particularly apparent when one carefully 
examines a standard ingredient in the reaction rate in use in the chemistry commu- 
nity - the flux-flux autocorrelation function for which we show that the use of our 
dividing surface and phase space approach allows the computation of this function 
without the long time integration of trajectories. 

In summary, our work on the geometry of reaction dynamics allows for a careful 
analysis and realization of Wigner's |Wig38] dynamical version for transition state 
theory. The dynamical foundations of Wigner's transition state theory did receive a 
great deal of attention in the 1970's in a series of seminal papers by Child, McLaf- 
ferty, Pechukas and Pollak [PM731 iPPTTl [PP78l TPTM IPP79bl lPCP80l lPC80l KTM, 
IPecBl) . and there is a wealth of dynamical ideas in these works. However, it is im- 
portant to realize that these works focus almost entirely on 2 DoF, and most of the 
results have not been generalized to 3 or more DoF. Nevertheless, for 2 DoFthey 
show how to construct a dividing surface without recrossing from the projection of a 
periodic orbit, the Lyapunov orbit associated with a saddle equilibrium point, to con- 
figuration space - the so called periodic orbit dividing surface (PODS) |PM73tlPP78j . 
In addition to this construction being limited to 2 DoFsystems, the Hamiltonian must 
be of type 'kinetic plus potential' - Coriolis terms due to a rotating coordinate system 
or a magnetic field are not allowed. 

The generalization to more than two degrees of freedom and to more general 
Hamiltonians has posed a major problem for decades. The reasons for the problems 
are twofold. On the one hand, a construction based on configuration space, as in the 
case of the PODS, simply does not work for systems with more than two degrees of 
freedom, as discussed in |WW04j . On the other hand, it was not clear what replaces 
the periodic orbit in higher dimension. For more than two degrees of freedom a 



periodic orbit lacks sufficient dimensionality to serve as a building block for the 
construction of a dividing surface. In fact, a completely new object, a so called 
normally hyperbolic invariant manifold (NHIM) |Wig94| takes the place of the periodic 
orbit in two degrees of freedom. It is interesting to recall a remark of Pechukas from 
his influential review paper |Pec81j : 



It is easy to guess that generalized transition states in problems with more 
degrees of freedom must be unstable invariant classical manifolds of the 
appropriate dimension, but to our knowledge no calculations have been 
done. 



Our work gives a precise characterization of these invariant manifolds in terms of 
the NHIM, as well as shows exactly what calculations are required to realize them 
in specific systems0 The NHIM is not only the building block for the construction 
of a dividing surface in arbitrary dimension, but it also forms the basis for locating 
the transition pathways for reactions in terms of the stable and unstable manifolds 
of the NHIM |WBW04bj . 

Finally, we remarked earlier that one way in which TST can "break down" is 
through deformation of the energy surface. Another way in which it may break down 
is through bifurcation of the NHIM. For 2 DoFsystems the NHIM is a periodic orbit 
and bifurcation theory for periodic orbits in Hamiltonian systems is well- developed 
( |Mey70[ IMH92j ). Bifurcation of the NHIM in 2 DoFsystems can lead to stable 
motions that "trap" trajectories in the transition region. This has been observed 
in |CP80l IMM97| . At present there exists no general bifurcation theory for NHIMs 
in systems with d DoF, d > 3, and this poses a limitation to the range of validity 
of our approach. In this case the relevant NHIMs are (2d — 3)-dimensional and 
contain their own nontrivial dynamics. The development of bifurcation theory for 
such objects promises to be a challenging and interesting mathematical problem that 
should yield new insights into reaction dynamics. 



The Dynamical Reaction Path: Thus far we have described the geometry of 
the energy surface near a saddle and the nature of the dividing surface that separates 
the energy surface near the saddle into two regions corresponding to reactants and 
products. Now we want to describe in more detail how trajectories approach, and 
move away from, the dividing surface. For this purpose the notion of the reaction 
path arises. 

Traditionally, the reaction path of a polyatomic molecule is the steepest descent 
path on the potential energy surface (if mass-weighted Cartesian coordinates are 

^It is perhaps worth pointing out that when reading the chemistry literature mathematicians might 
experience some confiision surrounding the phrases "transition state" and "dividing surface" . In some 
parts of the literature they are used synonymously. In other parts, they have a very different meaning, 
as can be seen from the above quote of Pechukas. A dividing surface cannot be an invariant manifold, or 
else trajectories could not cross the surface (trajectories on an invariant manifold remain on that manifold 
for all time). The confusion probably arose out of the PODS theory. In that situation the dividing 
surface and the invariant manifold (the periodic orbit) project to the same line in configuration space. 
The projection of a reactive trajectory to configuration space intersects this line in configuration space. In 
the three-dimensional energy surface however, the trajectory intersects the dividing surface and not the 
periodic orbit. The dividing surface is a 2-dimensional sphere, S"^, in this case (i.e. it is of one dimension 
less than the three-dimensional energy surface) and the periodic orbit is an invariant one-dimensional 
sphere, S^, that forms the equator of the sphere. The same situation holds for more than two DoF. The 
equator of our dividing surface is a normally hyperbolic invariant manifold (but a periodic orbit does not 
have sufficient dimensions to satisfy this requirement). 



used) connecting saddle points and minima ( [MHA80] ). Hence, it is a configuration 
space notion derived from properties of the potential energy surface that is used to 
describe a specific dynamical phenomenon. Similarly to TST, the literature related 
to reaction paths is vast. Searching ISI Web of Knowledge on the phrase "reaction 
path" yields more than 6,600 hits. Searching Google with the same phrase gives 
more than 7,900,000 hits. It is often assumed that a reacting trajectory, when pro- 
jected into configuration space, will be "close" to this reaction path, and much work is 
concerned with developing configuration space coordinates (and their associated con- 
jugate momenta) in which the dynamical equations that describe evolution "close" 
to this reaction path can be expr essed (see, e.g., |Mar66al [M^?66bl IMar68l IMHA801 
lMil83l INat92bl INat92al InGT+91[ iN^tOTl INat04l IHUBOTI IGGBOll IGBOSj ) . However, 
despite its fundamental importance in the historical development of the subject of 
reaction dynamics, one might question the relationship of this configuration based 
reaction path to the actual path taken in the course of the dynamical evolution from 
reactants to products. In fact, in recent years numerous experiments have shown that 
the actual dynamics may exhibit significant deviations from the "classical reaction 
path" |PCC+05l [SSHO2I IAYAD031 lLCZ+071 rTLL+04[ iBoOT IHK061 [PMnE 06] . 

In this paper we show that the coordinates given by normal form theory also 
give rise to an intrinsic dynamical reaction path, which is a trajectory on the energy 
surface. Its construction follows from the dynamical properties associated with the 
NHIM ("activated complex"). The NHIM has stable and unstable manifolds which 
as we will explain in detail have the structure of spherical cylinders, S'^'^~^ x R, 
and form the phase space conduits for reaction in the sense that they enclose the 
reactive trajectories. Our dynamical reaction path forms the centre line of these 
spherical cylinders, and gives rise to a phase space description of an invariant "modal 
partitioning" along the reaction path corresponding to energy in the reacting mode 
and energies in the (nonlinear) vibrational modes normal to the reaction path. 

1.2 Transition State Theory: Quantum Dynamics 

Historically a great deal of effort - mostly in the chemistry community - has been 
devoted to developing a quantum mechanical version of transition state theory (see 
the work by Miller and coworkers |Mil98aj ). Nonetheless, a quantum mechanical 
formulation of transition state theory is still considered an open problem (see the 
recent review by Pollak and Talkner |PT05b| ). The nature of the difficulties are 
summed up succinctly by Miller |Mil98b] : 

— the conclusion of it all is that there is no uniquely well-defined quantum 
version of TST in the sense that there is in classical mechanics. This is 
because tunnelling along the reaction coordinate necessarily requires one to 
solve the (quantum) dynamics for some finite region about the TS dividing 
surface, and if one does this fully quantum mechanically there is no 'theory' 
left, i.e., one has a full dimensional quantum treatment which is ipso facto 
exact, a quantum simulation. 

Part of the problem leading to this statement originates from (classical) transition 
state theory where the necessary theoretical framework to realise transition state 
theory for multi-dimensional systems as described in this paper has been developed 
only very recently. In particular, this realisation of classical TST requires one to 
work in phase space (as opposed to configuration space). This also has consequences 
for the development of a quantum version of transition state theory (which should 
reduce to classical TST in the classical limit and have all the computational benefits 



of a "local" theory like in the classical case). Again due to the lack of a theoretical 
framework, most approaches to developing a quantum version of transition state 
theory involve attempts to achieve a separation of the Schrodinger equation that 
describes the chemical reaction. However, like in the classical case this separation 
does not exist. In contrast to this, we will develop a quantum version of TST which 
is built in a systematic way on the classical theory presented in this paper. 
In the classical case the key idea to realise TST is to transform the Hamilton func- 
tion describing the reaction to normal form. In the quantum case we will establish a 
quantum version of the classical normal form theory, and from this all of the quan- 
tum reaction dynamics quantities will flow. In particular, the classical phase space 
structures that we found will play a central role in the computation of quantum 
mechanical reaction quantities. Quantum mechanical computations are notable for 
suffering from the "curse of dimensionality." We will see that the property of integra- 
bility which follows from the normal form in the classical case will have a quantum 
manifestation that renders computations of "local" reaction quantities tractable for 
high dimensional systems. This leads to very efficient algorithms for computing, e.g., 
quantum mechanical cumulative reaction probabilities and resonances. 

Classical normal form theory is a standard technique of dynamical systems the- 
ory, and there are many textbooks and tutorial articles that describe the subject. 
However, quantum normal form theory is probably much less familiar in both the 
dynamical systems community as well as the chemistry community. It is therefore 
useful to provide a discussion of the background, context, and historical development 
of the subject. 

Symplectic transformations like those involved in the classical normal form theory 
also have a long history in the study of Partial Differential Equations. In the theory 
of microlocal analysis they form one of the core techniques introduced in the late 60's 
and early 70's in the fundamental papers by Egorov, Hormander and Duistermaat, 
Ego69t IHor7H IDH72j . These ideas lead naturally to the consideration of normal 



forms for partial differential equations, and these were used to study the solvability 
and the singularities of solutions. The basic idea is the following. One can associate 
with a linear partial differential operator a function on phase space by substituting 
momenta for the partial differentials. The resulting function is called the symbol of 
the operator. One can now use a symplectic transformation to find coordinates in 
which the symbol has a particularly simple form. The crucial point now is that the 
tools from microlocal analysis allow one to quantise such a symplectic transformation. 
The result is a unitary operator which is called a Fourier integral operator and yields 
the transformation of the original partial differential operator corresponding to the 
symplectic transformation of its symbol (plus small error terms). This is the content 



of Egorov's Theorem, Ego69 . If the transformed symbol assumes a simple form, 
then the the transformed operator assumes a simple form too and its properties can 
be studied more easily. This construction was the basis for many developments in 
the theory of linear partial differential equations in the 70's, such as the study of 
the solvability and the propagation of singularities (see, e.g., the compendium by 
Hormander [Hor85al IHor85bj ) . 

In quantum mechanics the relation between operators and symbols mentioned 
above is the relation between the Hamilton operator, which defines a quantum me- 
chanical system, and the corresponding classical Hamilton function, which defines 
the classical dynamical system corresponding to the quantum system. The operator 
thus is the quantisation of the symbol, and microlocal analysis provides us with a 
powerful set of tools to analyse quantum systems. These ideas were applied, e.g., 
in the seminal work by Colin de Verdiere on modes and quasimodes |CdV77| where 



he constructed classical and quantum normal forms around invariant tori in phase 
space which still is an active area of research (see, e.g., the recent work by Cargo et 
al. [CGSL+n5p . 

In transition state theory the classical Hamiltonian relevant for reaction type dy- 
namics has an equilibrium point, and as we have discussed in the first part of the 
introduction one can use symplectic transformations to bring the Hamilton func- 
tion to a normal form in the neighbourhood of the equilibrium point. The tools 
from microlocal analysis will allow us to quantize this symplectic transformation 
and bring the Hamilton operator into a normal form, too. The problem of quan- 
tum normal forms near equilibrium points of the symbol has been studied quite 
extensively already. But most of this work concerns stable equilibrium points (see 
|BV901 IEGH9H |Sjo92[ IBGP99j ). Here the aim is to construct a quantum normal 
form in order to study energy spectra and eigenfunctions with very high precision. 
In the physics literature |Rob84l IAli85[ IEck86[ IFEBSal IFE88bj the same question 
was studied based on the Lie approach to classical normal forms. In the early works 
there has been some confusion about the ordering problem in quantisation, but these 
problems have been resolved by Crehan |Cre90j . 

The case of an unstable equilibrium point (or more precisely an equilibrium of 
saddle-center-- • - -center type), which occurs in transition state theory, has received 
much less attention in the literature so far. In this case one expects the operator to 
have continuous spectrum, and so instead of computing eigenvalues one is looking 
for resonances. Resonances are complex eigenvalues. Their imaginary parts are 
related to the finite lifetime of quantum states in the neighbourhood of the unstable 
equilibrium point. Since the problem is no longer selfadjoint, the determination of 
resonances is in general a much more difficult problem than that of eigenvalues (see 
|Zwo99j for a review). The case of a Hamilton operator for which the symbol has an 
unstable equilibrium is one of the few cases where resonances can be computed to 
high accuracy using a complex Bohr-Sommerfeld quantisation. For 2 DoFsystems, 
this was developed in |GS871 Sjo03) (for references in the chemistry literature, see, 
e.g., [SM91l[Moi98j where resonances are known as Gamov-Siegert eigenvalues). The 
methods were then extended to systems with more DoFby Sjostrand in |Sjo87| , and 
building on this work more complete results were obtained by Kaidi and Kerdelhue 
[KKOOj who derived quantisation conditions for the resonances which are valid to all 
orders in the semiclassical parameter h and are based on a quantum normal form. 
In |IS02] this was embedded into the study of more general normal forms for Fourier 
integral operators. 

The development and study of the quantum normal form near an equilibrium 
point of saddle-center- • • -center type is one of the mains aims in the quantum part 
of this paper. As mentioned above, the quantum normal form has already been used 
to study resonances in the literature before. We will see that the quantum normal 
form provides us with much more information which includes cumulative reaction 
probabilities and a detailed understanding of the dynamical mechanism of quantum 
reactions. To this end we will relate the quantum states described by the quantum 
normal form to the phase space structures that control classical reaction dynamics. 
In the classical case the NHIM is the manifestation of the activated complex. Due to 
the Heisenberg uncertainty relation, there is no such invariant structure in quantum 
mechanics. In fact the resonances will describe how the quantum activated complex 
decays. 

In order to use the quantum normal form to study concrete chemical reactions we 
have to be able to compute it explicitly, i.e., we need an explicit algorithm analogously 
to the classical normal form. The mathematical treatments in |Sjo87[ IKKOO] do 



not give us such an algorithm. Therefore we devefop a quantized version of the 
classical normal form algorithm which is similar to quantum normal form for stable 
equilibrium points in [CreOOl lEGHQlt IBGP99] . We give a complete exposition of 
our algorithm to compute the quantum normal form. At the level of symbols, the 
classical and quantum normal form algorithms are almost identical. The essential 
differences are that the Poisson bracket is replaced by the Moyal bracket, and rather 
than dealing with polynomial functions of the phase space coordinates, we deal with 
polynomial functions of the phase space coordinates and h. 

The outline of this paper is as follows. In Sec. [2] we will start by reviewing classical 
normal form theory. We will show in detail how to construct symplectic transforma- 
tions from the flows of Hamiltonian vector fields. The theory is presented in such 
a way that it allows for a direct comparison to the quantum normal form that we 
develop in Sec. [3l This section includes a careful review of the necessary tools from 
the symbol calculus which are required to quantize symplectic transformations. In 
Sec. m we discuss the phase space structures which govern classical reaction dynam- 
ics and show how these phase space structures can be realised with the help of the 
classical normal form. This includes the construction of a dividing surface, the role 
of the NHIM and its stable and unstable manifolds, the foliation of the NHIM by 
invariant tori and its relation to the activated complex, the definition of dynamical 
reaction paths and a formula for the directional flux through the dividing surface. 
In this section we will also relate the theory presented to the flux-flux autocorrela- 
tion function formalism that can be found in the chemistry literature. Sec. is the 
quantum mechanical analogue of Sec. [H We here use the quantum normal form to 
study quantum reaction dynamics. We show how to construct a local S-matrix from 
the quantum normal form and how this leads to an efficient algorithm to compute 
the cumulative reaction probability (the quantum analogue of the classical flux) . We 
study the distributions of the scattering states in phase space and relate them to 
the phase space structures governing classical reaction dynamics. We here will also 
relate the quantum normal form computation of the cumulative reaction probability 
to the quantum version of the flux-flux autocorrelation function formalism. In Sec. [6] 
we study quantum resonances that correspond to the (classical) activated complex. 
We will show how the resonances describe the quantum mechanical lifetimes of the 
activated complex. We will study the phase space distributions of the corresponding 
resonance states and interpret these distributions in terms of the phase space struc- 
tures associated with the classical dynamics of reactions. In Sec. [7] we illustrate the 
efficiency of the classical and quantum normal form algorithms for computing fluxes, 
cumulative reaction probabilities and resonances by applying the theory presented 
to several examples with one, two and three degrees of freedom. 



2 Classical Normal Form Theory 

In this section we summarise the main elements of classical Poincare-Birkhoff normal 
form theory for Hamiltonian functions. This is a well-known theory and has been 
the subject of many review papers and books [Dep69[ \DF76\ IAKN881 IMH921 IMur03j . 
The main reason for summarising the essential results here is so that the reader can 
clearly see the classical and quantum normal form theories "side-by-side" . In this way 
the classical-quantum correspondence is most apparent. This is explicitly illustrated 
by developing the classical normal form theory in a way that is rather different than 
that found in the literature. This difference allows us to explicitly show that the 
structure of the classical and the quantum normal form theories is very similar. At 



the same time, we emphasize that the classical normal form theory is an essential 
tool for both discovering and computing the necessary geometric structures in phase 
space with which we construct our phase space transition state theory in Section [H 
This section is organised as follows. In Sec. 12.11 we show how functions on phase 
space transform under symplectic coordinate transformations, which are constructed 
as Hamiltonian flows. In Sec. l2.2l we define what a (classical) normal form is and show 
how the formalism developed in Sec. 12.11 can be used to transform a Hamiltonian 
function into normal form to any desired order of its Taylor expansion about an 
equilibrium point. The general scheme is discussed in detail in Sec. 12.31 for the case 
of a saddle-centre-- • • -centre equilibrium point. 

2.1 Transformation of Phase Space Functions through 
Symplectic Coordinate Transformations 

The essence of classical normal form theory is to find a new set of coordinates, i.e., 
a change of variables, that transforms the Hamiltonian to a "simpler" form (and we 
will explicitly define what we mean by "simpler" shortly). Since we are dealing with 
Hamiltonian functions we want the coordinate transformation to preserve the Hamil- 
tonian structure, and this will be accomplished if the transformation is symplectic 
( |Arn781 IAM78j ) . A standard approach to constructing symplectic transformations 
is through the use of Lie transforms (see, e.g., |Mur03] ). which we now review. Before 
proceeding we note that there are issues related to differentiability of functions, exis- 
tence and uniqueness of solutions of ordinary differential equations, etc. However, we 
will proceed formally and assume that our functions have as many derivatives as re- 
quired and that solutions of ordinary differential equations exist, and are sufficiently 
differentiable, on domains of interest. Our purpose here is to develop methods and 
an algorithm. Its applicability must be verified for specific problems. 

A function W on phase space x M!^ defines a Hamiltonian vector field 



and at a point z = {q,p) = {qi, ■ ■ ■ jQcItPi, ■ ■ ■ iPd) in phase space this vector field 
takes the value 



The solutions of the ordinary differential equation ("Hamilton's equations") 



defines a Hamiltonian flow, z i— > z{e) := ^^r{z), which satisfies the properties 




(2.1) 




(2.2) 



(2.3) 




• ^"w = id, 



where id denotes the identity map, and 




(2.4) 



Most importantly for us, the Hamiltonian flow <I>^ defines a symplectic, or 'canon- 
ical', coordinate transformation of the phase space onto itself [Arn78j . This is signifi- 
cant because symplectic coordinate transformations preserve the Hamiltonian struc- 
ture. The Hamiltonian W is referred to as the generating function for the symplectic 
transformation 

We now consider the transformation of a (scalar valued) function on phase space 
under such a symplectic transformation. More precisely, for a phase space function 
A and a symplectic coordinate transformation defined from the flow generated by 
Hamilton's equations z{e) = ^^^(z), the transformation of the function under this 
symplectic transformation is given by 

A{e) = Ao^^-, (2.5) 

or, in coordinates, 

A{e){z{e))=A{z). (2.6) 

For our purposes we want to develop A{e) as a (formal) power series in e. We 
begin by computing the first derivative of A{e) with respect to e giving 

^^(e) = -{VA, Xw) o = {W, A} o (D^^ , (2.7) 

where VA = (^§^, . . . , . . . , ^) is the gradient of A, (•, •) is the standard 

scalar product in R'^'^, and 

is the Poisson bracket of W and A. Using the fact that W is invariant under the flow 
( i.e., W {^^) = W (<^vf)' ^'-i ™ other words, the Hamiltonian W is constant 
along trajectories of the vector field X\y generated by W) we can rewrite (12. TP as 

^A{e) = {W,A{e)}. (2.9) 

The Poisson bracket gives us a convenient way of representing the derivatives of a 
function along trajectories of Hamilton's equations. We simplify the notation further 
by defining the adjoint operator 

adiy : A ^ adw A := {W, A} (2.10) 



associated with a generating function W. We can now differentiate ()2.9p again to 
obtain the second order derivative with respect to e. 



de2 



^^'^ = ^(^^^'0 =W^^(^)} = i^'{^'^(^)»=^ [^dw]'A{e). (2.11) 
Continuing this procedure for higher order derivatives gives 



d" d / d"~i \ d 

^nMe)-^^{^Aie))^{WA-{W,-Aie)}-}} 



-.{W,{---{W,{W,Aie)}}---}} 
--: [adi^]"^(e). 



(2.12) 



Using these results, we obtain the Taylor expansion of ^(e) about e = 0, 



^(^) = E^d^^(^)Uo = E^[-dH.]"^, (2.13) 

n=0 ■ n=0 

where ^(0) = A and [ adw ] are defined as in Equations ()2.9p - ()2.12p with [ adw ] 
A. 

Equation ()2.13p gives the Taylor expansion with respect to the flow parameter or 
'time' e for a phase space function A that is transformed by a symplectic transforma- 
tion defined by the Hamiltonian flow generated by the function W. It will form the 
basis of the classical normal form method where the idea is to "simplify" (or "nor- 
malise") a function which, for us, will be a specific Hamiltonian through the choice of 
an "appropriately chosen" sequence of symplectic transformations that simplify the 
Hamiltonian "order by order" of its Taylor expansion with respect to the phase space 
coordinates z = {q,p). First, we need to make clear that the normal form procedure 
that we develop here is valid in a neighborhood of an equilibrium point. This means 
that the normal form is a local object whose dynamics have meaning for the original 
Hamiltonian only in a neighborhood of an equilibrium point. In order to describe 
the terms in the Taylor expansion of a given order in the phase space coordinates 
more precisely we introduce the vector spaces W^j, s G Mq, of polynomials which are 
homogeneous of order s. The space W^j is spanned (over C) by all monomials of the 
form 

d d 
q»pP := Yl ql'^pl" , where \a\ + \(3\ ■.= Y,o^k + (ik = s . (2.14) 

k=l k=l 

The following two lemmata are the key tools used in the computation of the 
classical normal form. 

Lemma 1. LetW £ W^[, A £ W^i with s,s' >1, then 

{W,A} (2.15) 

and for n > 0, 

[adH^]"^G (2.16) 
if n{s' — 2) + s > and [advi/]"^ = otherwise. 

Proof. This Lemma can be proven by direct calculation. □ 

This lemma is key to the proof of 
Lemma 2. Let W £ with s' > 3 and 



A = Y,^s (2.17) 



with As G W^V Then 



oo ^ oo 

A' ■.= Ao<^>^^ = Y^-[adwrA = Y,K , (2.18) 



n=0 s=0 



where 



A's= ;^[adiy]"^-n(.'-2), (2.19) 

n=0 

where [-p^] denotes the integer part of . 



Proof. Using (12.170 . we write out the next to last term in (12.180 as a series of series 
as follows (where we have also changed the summation index from s to j in order to 
avoid possible confusion): 



n=0 n=0 j=0 

oo oo oo ^ 



2 

j=0 j=0 j=0 

oo ^ oo ^ 

+ E^[adH^]% + --- + E;^[ad;^]% + --- •(2-20) 

j=0 ■ j=0 

We now want to inspect each series in the series and extract the order s term from 
each one. Then summing these terms will give the series (j2.19p . Using Lemma [H we 
find 

[adH/]%eW^'-')+^ (2.21) 
Now we wish to choose j such that 

[adw^Aj eW!i. (2.22) 
Comparing (12.210 and ()2.22p . this is true for 

j = s-n{s'-2). (2.23) 

Hence it follows that 



A's= E ^[^dwrA-nis'-2). (2.24) 

□ 



n! 

n=0 



2.2 Definition and Computation of the Classical Normal 
Form 

We will now define when a Hamilton function is in classical normal form. Here 
we use the adjective 'classical' to distinguish the normal form in the case of classical 
mechanics from the normal form that we will define for the case of quantum mechanics 
in Sec. El As we will see, in general a Hamilton function is not in normal form. 
However, as we will show in detail, the formalism reviewed in the previous section 
can be used to construct an explicit algorithm which allows one to transform a 
Hamilton function to normal form to any desired order of its Taylor expansion. 

The starting point is a Hamilton function with an equilibrium point at z = zq, 
i.e., VH{zo) = 0. Let H2{z) := ^{z — zq,!)'^ H{zo){z — zq)) be the quadratic part of 
the Taylor expansion of H about zqH We then make the following 

Definition 1. We say that H is in classical normal form with respect to the 
equilibrium point zq if 

adn.H = {H2,H} = 0. (2.25) 
^Here , D^iJ(zo) denotes the Hessian of H at zq, i.e., the matrix of second derivatives {dzidzjH{zo))i 



It follows from this definition that if H is in normal form then H2 will be an 
integral of the motion generated by the Hamilton function H and moreover, as we 
will see below, depending on the structure of H2, further integrals of the motion 
will exist. A consequence of the existence of integrals of motion is the structuring, 
or foliation, of the phase space by lower dimensional surfaces or manifolds that are 
invariant under the dynamics. If we choose initial conditions for Hamilton's equations 
then these initial conditions will determine values of the integrals of motion. The 
full solution of Hamilton's equation will then be contained in the manifold given by 
the common level set of the integrals corresponding to the initial values. This way 
the integrals of the motion confine the possible dynamics. Moreover, the existence 
of integrals of the motion significantly simplifies the study of the dynamics. 

In general a Hamilton function is not in normal form. However, we will use the 
formalism and results developed in the previous section to transform a Hamiltonian 
to normal form in a neighbourhood of the equilibrium point to a certain order of 
its Taylor expansion about the equilibrium point. As we will see, the transformed 
Hamiltonian function truncated at this order will lead to a very accurate description 
of the motion in the neighbourhood of the equilibrium point. What we mean by 
"accurate description" is discussed in Sec l4.5p . 

We develop the following procedure. We begin with our "original Hamiltonian" 

H = , (2.26) 

and we construct a consecutive sequence of symplectic transformations 

i^(0) ^ ^(1) ^ ^(2) ^ ^(3) ^ . . . ^ H^N) ^ (2.27) 

where is a sufficiently large integer which will be the order at which we will truncate 
the normal form series. 

The first step in the sequence (|2.27p is obtained by shifting the critical point zq 
to the origin of a new coordinate system. We set 

^(1) =z-zo . (2.28) 

The Hamiltonian function H^^^ is the representation of H'^^'^ in terms of the new 
coordinates z^^\ i.e., 

i?«(z«) = i/(o)(z«+zo). (2.29) 

Once the equilibrium point is shifted to the origin, our normal form procedure 
will require us to work with the Taylor expansion of the Hamiltonian H^^^ about the 
origin in a "term-by-term" fashion. Let 

00 

H^^^ =Eo + Y,H'^^\ (2-30) 

s=2 

where 

H^^Kq.p):= -^d^d^H('H0,0)q''p^ (2.31) 

\a\ + \l3\=s 

are the terms of order s. Here we employ the usual multi-index notation; for a = 
(ai, . . . , Ofi) € Nq we have \a\ = ai + . . . + a^, al = ai!a2' • • • 9° = Qi^Q2^ ■ ■ ■ Q^'^ 
and dg = -^n^ ■ ■ ■ (for /3 G Nq and p G W^, the notation is analogous). Since 

()2.30p is a Taylor expansion of a Hamiltonian about an equilibrium point at the origin 



it follows that H\ =0. In particular, Hq = Eq is the "energy" of the equilibrium 
point. 

At the next step in the sequence ()2.27p we choose a linear symplectic transfor- 

(2) 

mation such that H2 assumes a "simple form". In other words, we seek a trans- 
formation that simplifies the quadratic part of the Hamiltonian or, equivalently, the 
linear part of the Hamiltonian vector field. This is accomplished by choosing an 
appropriate symplectic 2d x 2d matrix M, i.e., a matrix statisfying M'^ J M = J, 
where J is the standard 2d x 2d symplectic matrix 

whose blocks consist of d x d zero matrices and d x d identity matrices. We then set 

z(2) = MzW, (2.33) 
and the corresponding transformed Hamiltonian is given by 

^{2)(^{2))^^{l)(^^-l^(2))^ (2.34) 

(2) 

Which form of H2 can be considered to be "simple" depends on the nature 
of the particular equilibrium point (i.e., the eigenvalues and eigenvectors associated 

with the matrix obtained by linearising Hamilton's equations about the origin). The 

(2) 

main benefit of having H2 in a "simple" form is that this will simplify the explicit 
implementation of the algorithm to normalise the higher order terms, n > 3, i.e., 
how to choose the next steps in the sequence ()2.27p . Therefore, "simplify" could 
mean that we would seek a transformation that would diagonalise the linear part of 
Hamilton's equations, or transform it to "real Jordan canonical form" in the case 
of complex eigenvalues. Clearly, constructing such a transformation is a problem 
in linear algebra for which there is a large literature. However, the symplectic case 
tends to bring with it new difficulties, both in the analytical and computational areas 
(see, e.g., |CK99] ). In the next section we will see how to simplify the linear part 
of Hamilton's equations for our particular case of interest, i.e., a saddle-centre- • •- 
centre equilibrium point satisfying a certain "nonresonance" condition. However, it is 
important to realise that the normal form algorithm does not depend on the specific 
form taken by the linear part of Hamilton's equations. 

Up to this point we have located an equilibrium point of interest, translated it 
to the origin, Taylor expanded the resulting transformed Hamiltonian H^^^ about 
the origin (for which h[^^ = 0), and constructed a linear symplectic transformation 

in such a way that the quadratic part of the resulting transformed Hamiltonian, 

(2) 

Hp, is "simple". Now we are ready to describe how to normalise the terms of 
order three and higher, i.e., how to define the next steps in the sequence ()2.27p . To 
accomplish these transformations we will use the formalism reviewed in Sec. 12.11 and 
successively transform the Hamiltonian by the time one maps of the flows generated 
by Hamiltonian vector fields. More precisely, for n > 3, H^"'^ is computed from 
according to 

00 ^ 

H(n) ^ ^(n^l) o CD^^ = _ [ad^j'^i^(n-l) (2.35) 

k=0 

with a generating function Wn G W^. The order s term of the Taylor expansion 
of H^"'^ expressed as a series involving terms in the Taylor expansion of //("-I) and 



Wn is obtained by substituting the Taylor expansion of into ((2:35]) and usme 

Lemma [5J This gives 

^i"^ = E fei [ ^d^" ] 'Ht~Hl2) > « > 3 . (2.36) 

/c=0 

The corresponding transformation of phase space coordinates is then given by 

= $il,Jz("-i)), n>3. (2.37) 

We note that in fact also the afiine linear symplectic coordinate transformations 
()2.28p and (j2.33p which formed the first two steps in the sequence ()2.27p can be 
formally expressed as time one maps of Hamiltonian flows with generating functions 
Wi G and W2 S W^j, respectively. A generating functions Wi whose time one 
map achieves the translation ()2.28p is given by 

Wi{z) = -{zo,Jz), (2.38) 

where J is the standard 2d x 2d symplectic matrix defined in Equation ()2.32p . This 
gives 

^L{z)=z-zo. (2.39) 



z 



(1) 



In this case the upper limit of the sum in ()2.36p is infinity. It is in general not 
straightforward to explicitly give an expression for a generating function W2 G W^j 
whose time one map achieves the linear symplectic transformation (I2.33P for a given 
symplectic matrix M. But such a W2 always exist^. For n = 2 in Equation (j2.36p the 
upper limit of the sum is again infinity. In the next section we will provide a matrix 
M which achieves the simplification of the quadratic part of the Hamiltonian function 
for the case of a saddle-centre-- • • -centre equilibrium point satisfying a nonresonance 
condition without specifying the corresponding W2- Note however that it is M and 
not necessarily W2 which is required for our normalisation procedure. 

Let us now proceed with the nonlinear symplectic transformations generated by 
polynomials Wn S W"[ with n > 3 to achieve the third and higher steps in the 
sequence ()2.27p . The first thing to note is that these transformations will not alter 
the zeroth order term, £^0) and we will also have h[^^ = h["'^ = 0, n > 3. The zeroth 
order term is unaltered since the upper limit in the sum ()2.36p is zero for s = 0. The 
first order term stays zero because for s < 1 in combination with n > 3 and s = in 

combination with n = 3, the upper limit in the sum ()2.36p is again zero. For n = 3 

(2) 

in combination with s = 1, the upper limit is 1. However, the k = 1 term, adw-j, Hq , 
in the sum (j2.36p is zero because Hq is the constant Eq and hence vanishes when 
adws is applied to it. 

(2) 

Moreover, the quadratic part of the Hamiltonian H2 will not be modified by 
the transformations generated by Wn, n > 3. We will show this directly from our 
formalism. 

Lemma 3. ijj"^ = H^^\ n > 3. 



^This follows from two facts. Firstly, the group of linear symplectic transformations is connected, and 
therefore the image of the exponentiation of its Lie algebra is connected, too. Secondly this Lie algebra 
is isomorphic to the vector space of quadratic polynomials endowed with the Poisson bracket |Fol89| . 
Therefore the set of all time one maps generated by quadratic elements of W^j is the whole symplectic 
group. 



Proof. The idea is to use ()2.35p to transform from H^" to H^"'\ and then to show 
that //^"^ = //^""^^ for n > 3. 

We separate out the constant and quadratic parts of H^^~^^ as 

oo 

//("-I) =E^ + i?^"-') + (2.40) 
and then we substitute this into (|2.35|) to obtain 

oo ^ oo ^ oo ^ oo 

A;=0 ■ A;=0 " A;=0 ' s=3 

(2.41) 

Note that the first series in this expression only admits the A; = term, Eq. We 
consider the case n > 3. In this case, the third series, using Lemma [U only admits 
terms of order larger than or equal to three. Hence, all of the quadratic terms must 
be in the second series. Using Lemma [H the /c*^ term in that series is contained in 
yyfc(?i 2)+2^ Therefore the only quadratic term occurs for A: = 0, which is ffg" 

□ 

Lemma [3] motivates the definition of the operator 

P:=ad^(2) ={/7f ,•}. (2.42) 

In fact, T) will simply be a convenient shorthand notation for the operator ad/^j = 
{i?2, •} in the definition of the the normal form in Definition [T] in terms of the coor- 
dinates z^'^\ The operator 2? plays a crucial role in the computation of the normal 
form transformation. 

The other important point to realise when transforming to H^'^^ with , 

Wn G VV^, is that all terms of order smaller than n are unchanged (however, the terms 
of order larger than n are modified by the n^^ order normalisation transformation). 
This is essential for the success of the iterative process and we provide a proof of this 
result now. 

Lemma 4. For n > 3 and < s < n, = . 

Proof. First, it is important to consider the upper limit of the sum ()2.36p . For 
0<s<n — 3itis zero, which indicates that for these values of s only the A; = 
term is nonzero. Hence, we have 

= , < s < n - 3. (2.43) 

Next we separately consider the cases s = n — 2 and s = n — 1. Using (I2.36P we find 
for s = n — 2, 

H^:\ = Ht~2^ + adH^„ Ht'^ = Ht,'^ (2.44) 
since h'^ = Eq = const.. For s = n — 1, (|2.36p gives 

= Ht,'^ + adw„ Ht'^ + 5„,sl[^dw„]'Ht'^ = Ht\'^ (2.45) 

since = and H^^ = Eq = const.. The Kronecker symbol in the last term 

of the second expression shows that this term only occurs for n = 3. □ 



Now if we consider the rfi^ order term in this will show us how to choose 
Wn, n > 3. 

Lemma 5 (Homological Equation). For s = n > 3, 

i/W = Hi^n-l) _'pWn . (2.46) 

Proof. This result is also obtained from ()2.36p . with a careful consideration of the 
upper limit of the sum. The case n > 5 is the most straightforward. In this case only 
A; = and k = 1 contribute in the sum, and using (12. 8p . we obtain immediately that 

H^^^ = //("-I) + adv^„ ) = - ad^(.) Wn = //(""^^ -VWn- (2.47) 

The special cases s = n = 4 and s = n = 3 must be considered. These will give rise 
to some additional terms in (j2.36p . However, as for Lemma HI these will be zero if 
we take into account Hi =0 and [advK„ \ Eq = for integers A; > 0, n > 3. 

□ 

Equation ()2.47p is known as the homological equation. We want to solve the 
homological equation, i.e., find a function Wn E W^, in such a way that H^"^^ is in 
normal form up to order n. To this end note that it follows from Lemma [T] that T> 
defines a linear map of into W^, i.e., for each n, 

V : >V"i ^ >V"i. (2.48) 

In order to have H^"^^ in normal form up to order n we have to require T> = 0. 
Looking at the homological equation (j2.47p this means we need to find a function 
Wn e W"i such that H^^ = Ht~^^ - P VF„ is in the kernel of the restriction of V to 
W3, i.e., 

H^^^ =H^^-^^ -VWne KerPj (2.49) 

cl 

Definition 2. We will call the homological equation (|2.46p solvable if for any n>3 

there exist for any Hn G an Wn G W"i such that 

Hn-VWneKerV\^^^ . (2.50) 

Whether the Homological equation is solvable and how such a Wn can be found 
depends on the structure of P, i.e., on the structure of the matrix associated with the 
linearisation of the vector field about the equilibrium point. In the next subsection 
we will show that the homological equation is solvable in the case of a saddle-centre- 
• • • -centre equilibrium point and explain how Wn can be found. 

We summarise the results of this section in the following 

Theorem 1. Assume that a Hamiltonian function H has an equilibrium point at 
zq G R*^ X M,'^, and that the homological equation is solvable. Then for every N ^ 
there is a symplectic transformation <I>7v such that 

//ocD^i = /7(5Jp + 0;v+i, (2.51) 

where H^Qjqp is in normal form (with respect to z = (0,0) j and Otv+i is of order 
+ 1, i.e., there exists an open neighbourhood U of z = (0,0) and a constant c > 
such that 

|0^+i(6z)| <ce^+^ (2.52) 

for z G U and e < 1. 



Proof. Following the scheme described in this section we normalise the Hamilton 
function H order by order according to the sequence (j2.27p . We start by choosing a 
new coordinate system z^^^ = z — zq which has the equilibrium point zq at the origin 
(see (|2.28p ). and Taylor expand the Hamilton function H^^\ which we obtain from 
expressing H in the new coordinates z^^^ (see Equation ()2.29l) ). about z^^^ = to 
order N. The remainder which we denote by -R^v+i then of order + 1. 

We then choose a symplectic 2d x 2d matrix M to define a linear symplectic 
transformation to new coordinates z^^) = Mz^^^ in terms of which the quadratic 
part of the transformed Hamilton function H^"^^ (see Equation (|2.34|) ) assumes a 
simple form. As mentioned above the choice of M depends on the nature of the 
equilibrium point and will simplify the calculation of the steps for n > 3 in the 
sequence (I2.27p . Apart from this however, the choice of the symplectic matrix M is 
not important. We thus get 

N 

h('^=Eo + Y.hP+R^^I„ (2.53) 

s=2 

where hP{z^'^^) = hP{M-^ z^'^'>), i.e. hP £ W^^ ior s = 2, . . . , N, and the remain- 
der term R^^\i given by i?j^^-^(z(^)) = R^^^^^{M^^ z^^)) is again of order + 1. 

Having simplified the quadratic part, we proceed inductively by subsequently 
choosing generating functions Wn S VV^, which at each order n, n = 3, . . . , N , solve 
the homological equation (|2.46p and determing H^"'^ from as follows. For 

n > 3, H^"'^^^ is of the form 

N 

(2.54) 

where i/i" G W^j and R^^_^l^ is of order A'" + l. Using this decomposition of 
we can write for H^"-^ = o <I>^^, 

N 

^(n) ^ ^i,(-l)ocI,^l +i?(;;,l)o$^^ (2.55) 
s=0 

N oo 

s=0 k=0 
s=0 k=Q 

N oo ^ 

<i = <;i'^°^^l + E E y[-<^w.fH(--'K (2.58) 

We here used Eq. ()2.35p to get ()2.56p . To obtain ()2.57p from ()2.56p we removed all 
those terms from the double sum in ()2.56p contained in the W^j with s > + 1 

and absorbed them in the new remainder term p.58p . Since the symplectic 

transformations ^ly^^ are near identity transformations for n > 3 the remainder term 

-^TV+i again of order + 1. 

After the step n = N the terms of order less than or equal to A^ of the Hamilton 
function H^^^ are then in normal form (with respect to z = (0,0)). The symplectic 



where 



transformation <I>7v in Eq. p.5ip and the corresponding new coordinates z^^^ are 
then given by 

z(^) = $^(z) = cDi^^o...o$i^^(z(2)), z(2)=Mz«, z«=z-zo. (2.59) 

□ 

From the point of view of appUcations the definition of the normal form in Defini- 
tion [1] is not very practical since it requires one to carry out the procedure described in 
the proof of Theorem[T]for — > oo. In general, it is well known that such normal form 
transformations do not converge, except in special cases |SM71llBru71t[Rus67tlPM03] . 
For applications it is more practical to consider the truncated normal form. 

Definition 3 (A^**^ Order Classical Normal Form). Consider a Hamilton function H 
with an equilibrium point zq G M!^ x R"' which, for G IN, we normalise as described 
in TheoremUi Then we refer to H^^p in Equation (I2.5ip as the N^^ order classical 
normal form (CNF) of H. 

Note that in order to compute the A^**^ order normal form it is sufficient to carry 
out the Taylor expansion of the Hamiltonian up to order A^. The remainder term can 
be neglected immediately since the procedure described in the proof of Theorem [T] 
shows that no terms from the remainder term will enter the A^*'^ order normal form. 

Of course, the normal form procedure presented in this section raises questions 
like "what is the error associated with truncating the normal form at some finite 
order?" After all, one is interested in the dynamics associated with the full, origi- 
nal Hamiltonian. Another obvious question is "what is the optimum order at which 
to truncate the normal form so that errors are minimised?" There is no general 
theory that can be used to answer such questions. They must be addressed on a 
problem- by-problem basis. Fortunately, truncating the normal form does give ex- 
tremely accurate results in a number of problems |WBW04at IWBW04bl IWBWOSbj . 
and we will consider this in more detail in Section f4.5[ 

2.3 Nature and Computation of the Normal Form in a 
Neighborhood of an EquiUbrium Point of Saddle-Centre- 
• • • -Centre Stabihty Type 

We now describe the computation of the normal form in the classical situation of 
interest to us; in the neighborhood of an equilibrium point of saddle-centre-- •• - 
centre stability type. This means that the matrix associated with the linearisation 
of Hamilton's equations about the equilibrium point has two real eigenvalues, it A, 
and d—1 complex conjugate pairs of pure imaginary eigenvalues, ±icok, k = 2, . . . , d. 
Moreover, we will assume that the w^, k = 2, . . . ,d, are nonresonant in the sense 
that they are linearly independent over the integers, i.e., A;2 a;2 + • • • + ^cZ i^cZ 7^ for 
ah {k2,...,kd) £ - {0} (note that the more stringent diophantine condition for 
nonresonance ( [AKN88] ) is not required for our work). 

But first, we locate the equilibrium point of interest, denote it by zq = {qo,Po), 
and translate it to the origin using the generating function given in (j2.38p . The 
Taylor series of the corresponding Hamiltonian then has the form 

CO 

/7(i)(za)) = Eo + f(^)(z«) + ^ fW(z«). (2.60) 

s=3 



We next construct a linear symplectic transformation M : ]R^°' i-^ such that 
for z*^^) = M z^^\ we have 



Hf\z^-)) = Xpf\f^ + E Y(rf^)^ + if'?) ■ (2-61) 

k=2 

We note that for some purposes it is convenient to consider also a slightly modified 
version of the coordinates z^'^'^ = {q^'^\p^'^^) which for later reference we will we 
denote by (Q(2)^p{2))_ The coordinates (q^^^^j^^)) ^nd (Q(2)^p(2)) 

agree in the centre 

(2) (2) (2) (2) 

components, i.e., Q\ = q\ and — Pk ^or k = 2, . . . ,d, but are rotated versus 
each other by an angle of 45° in the saddle plane, i.e., 

Q? = ^^{<l?-P?), P?' = ^^{^?+P?). (2.62) 

Note that the tranformation from {q^'^\p^'^^) to {Q^'^\ P^'^^) is symplectic. In terms 
of {Q^'^\ P^'^^) the quadratic part of the Hamiltonian assumes the form 

i/f (Q(^P(^)) = ^((Pf))^ -Qf^)^) +X;T((^f + • (2-63) 

k=2 

The quadratic part then consists of the sum of one inverted harmononic oscillator 
(or "parabolic barrier") and d—1 harmonic oscillators. 

In order to construct the 2d x 2d matrix M above we label the eigenvalues of 
JT)^H{zo) (which is the matrix corresponding to the linearisation of Hamilton's 
vector field around the equilibrium point) in such a way that 

ei = -ei+d = A, efc = -efc+d = iwfc , k = 2,...,d, (2.64) 

and then use the corresponding eigenvectors vi, . . . ,V2d to form the columns of the 
matrix M according to 

M = {ciVi,C2Rev2, ■ ■ ■ ,CdRevd,ciVi+d,C2lmv2, ■ ■ ■ ,CdlniVd) , (2.65) 

where ci , . . . , q are scalars defined as 

:= (ui, Jwi+d) , c^"^ := {Revk,Jlmvk) , k = 2,...,d. (2.66) 

The constants ci, . . . , q guarantee that the matrix M will be symplectic, i.e., M will 
satisfy M^JM = J. Here we have assumed that the eigenvectors vi and vi^d have 
been chosen in such a way that {vi,Jvi+d) is positive (if {vi,Jvi+d) < then we 
multiply vi^d by -1). It is not difficult to see that c^^, k = 2, . . . ,d, are automatically 
positive if the frequencies uJk are positive 0. Using the fact that {vn, Jvk) = for n 
and k from the distinct sets {1,1 + d}, {2,2 + d}, {d, 2d}, it is easily verified 
that the matrix M satisfies M'^ J M = J. 



^ In fact, if one of the is negative then this means that the corresponding frequency is negative; this 
is a case which we have excluded, although it can be dealt with in a way that is similar to the procedure 
described in this paper. 



2.3.1 Solution of the homological equation 



Given a Hamiltonian function whose quadratic part is of the form ()2.6ip , the solution 
of the homological equation derived in Lemma [5] for any order n > 3 is extremely 
simple and transparent if we first perform the following symplectic complex linear 
change of coordinates z^"^ = {q^"'\p'^"'^) ^ {x, ^) which has the components xi = 
^1 = and 



-k:=^{cit^-wt^), ik:=^{pt^-nt^), fc = 2,...,d. (2.67) 

Here, and for the rest of this section, we omit the superscript (n) for x and ^ for the 
sake of a simpler and less cumbersome notation. 

In terms of the phase space coordinates (x,^), the linear map T) takes the form 

d 

V = A(^i% - xid-j,^) + ^icjfc(^fe% - Xkdxk) ■ (2.68) 

k=2 

The form of ()2.68p is significant for two reasons. One is that when the monomials of 
order n defined in (j2.14p are expressed in terms of the coordinates (x, ^) they form a 
basis for W^. We have 



= span I x'^e ■= J] xfit : |a| + := «fc + A = n i 

k=l k=l ^ 



(2.69) 



Secondly, in this basis the linear map (I2.68|) is diagonal. In fact, using ()2.68p . we see 
that the image under P of a monomial x"^^ G is 



u yr a \ d 

^\{<''€'' = [m-o^i) + Y.''Mh-o^k)) Wx^il". (2.70) 

fc=l ^ k=2 ' k=l 

These monomials thus are eigenvectors of (I2.68p . 

Since the map P can be diagonalised it follows in a trivial way that can be 
represented as the direct sum of the kernel of V acting on WIS, Ker PL.,„, and the 

image of V acting on W^, Im ^^|yyn! i-e., 

= Ker V\ © Im V\ (2.71) 

cl cl 

Now we can express i?n" as 

Ht-'^ = <Ker^ + (2.72) 
(n~l) I fri — 1) I 

where -ff^-xer ^ -^^^ -^Iw^ ^n-im ^ then choose Wn such 

' cl ' cl 

that 

VWn = hI^-J^ , (2.73) 

and therefore by ()2.46p 

The choice of Wn is not unique since one can always add terms from the kernel of 
P|^„. However, we will require Wn G Im I?|^„, i.e., we will invert V on its image 

cl cl 

Im 2^1^71, which renders the choice of Wn unique. 



Using our assumption that the frequencies iJ2,---,uJd are nonresonant, i.e., hn- 
early independent over we see from (|2.7U|) that a monomial x°'^^ is mapped to 
zero if and only if ak = (5k for all A; = 1, . . . , d. In particular Ker ^^|yys = {0} if s 

cl 

is odd. This implies that coordinate transformations can be constructed such that 
all odd order terms are eliminated. Moreover, for s even, the terms that cannot be 
eliminated are those which are sums of monomials for which Xj. and have equal 
integer exponents for all /c = 1, . . . , d. 

Concretely, we can compute Wn according to (j2.73p as follows. We assume that 
^n-lra linear combination of L monomials of order n. 



<im =E^'n-r€^^ (2-75) 

1=1 k=l 



with Ylk=i^k;i + Pk;i = u ioi sll I = 1, . . . , L, and for all / = there is 

at least one k = 1, . . . , d for which a^-i 7^ Pk;i (i-e-, the vectors {ai-i, . . . , Ud-i) and 
{Pi-i, . . . , /3d;/) are different for all / = 1, . . . , L). Upon inspecting ()2.70p . and using 
()2.73p . we see that a generating function Wn that solves the homological equation is 
given by 

wn=t vr, — . — T n <r €r • (2-76) 

1=1 - + Z^k=2 ^^k{Pk;l - ak-i) 

As mentioned above this solution of the homological equation is unique if we require 
Wn to be in Im ^?|yyn. 

2.3.2 Integrals of the classical motion from the A^**^ order classical 
normal form 

The A^*'^ order classical normal form is a polynomial in the A^**^ order phase 

space coordinates (I2.59P which, in order to keep the notation in this section sim- 
ple, we will denote by {q,p), i.e., we will omit superscripts (A^) on the phase space 
coordinates. As discussed in the previous section, it follows that if we perform the 
symplectic complex linear change of coordinates xi = qi, = pi and 

Xk-=^iqk-m), Ck ■■= -^{pk - m) , k = 2,...,d, (2.77) 

then the coordinate pairs Xk and will have equal integer exponents for all k = 
1, . . . , d in each monomial of H^-p. As a consequence the functions 

I = PiQi = (,ixi , J,^ = ^[pl + ql) =i^^.Xk, k = 2,...,d, (2.78) 

are integrals of the motion generated by H^p. This assertion is simple to verify 
with the following computations: 

^/ = {/,i7g;i,} = 0, Aj, = {jfc,i/g;i,} = 0, k = 2,...,d. (2.79) 



The integrals of the motion / and can be used to define action angle variables. 
We therefore define the conjugate angles 



^ I (2l±|i) , p,q, < 

I tanh-i (21^) , pigi > ' (2.80) 

= arg(pfc + igfc) , k = 2,...,d. 

It is not difficult to see that the map {q,p) ^ {(pi, . . . , fd, I,J2,---, Jd) is symplectic. 

For k = 2, . . . ,d, the ranges of the are [0, 2tt) and the ranges of the Jfc are 
[0, oo). The maps {qk,Pk) ^ ivk, Jk) are singular at = Pfc = where the angles ipk 
are not defined. Away from the singularities the maps are one to one. In contrast, 
the range of both ipi and / is R {ipi thus is not an angle in the usual sense). The 
map {qi^pi) [^i^ h) is singular on the lines pi = and qi = which map to I = 
with ifi = oo and (fi = — oo, respectively. Even away from the singularities each 
{ipi,I) has still two preimages {qi,Pi) which correspond to the two branches of the 
hyperbola / = piqi- The coordinate lines of the action angle variables are shown in 

Fig.m 

We note that in terms of the coordinates {Q,P) with {Qk,Pk) = {qkiPk)-, k = 
2, . . . ,d, and 

Ql = ^{qi-Pl), Pi = -^(gi+pi), (2.81) 
the integrals Jk, k = 2, . . . ,d, are of the same form while / changes to 

I = l{P?-Ql)- (2.82) 

The angles cpk, k = 1, . . . ,d, are cyclic, i.e., the Hamilton function -ff^NF ^ffec- 



later reference we introduce the function K^-p defined via 



tively depends only on the integrals / and Jk, k = 2, . . . ,d. To indicate this and for 
;nce we intro 

-"CNF ~ -^CNFV-' ' "^2, • • • , JdJ (2 83) 

= Eq + XI + u;2J2 + ■ ■ ■ + ^dJd + higher order terms . 

Here the higher order terms are of order greater than 1 and less than or equal to 
[A^/2] in the integrals, where [A^/2] denotes the integer part of N/2. Note that since 
the Hamilonian in normal form does not have any odd order terms, only the case of 
even A'^ is of interest. 

As we will see, the classical integrals of motion are extremely useful for character- 
ising, and realising, classical phase space structures. However, the obvious question 
arises, and must be answered. These are constants of the motion for the A^**^ order 
classical normal form H^p. How close to being constant are they on trajectories of 
the full Hamiltonian? Also, we will use them to construct certain invariant manifolds 
for the N^^ order classical normal form Hq^p- How close to being invariant will these 
manifolds be for the full Hamiltonian? These questions must be asked, and answered, 
on a problem- by-problem basis. A number of studies have recently shown that for 
moderate A^ (e.g. 10-14), these integrals are "very close" to constant for the full 
Hamiltonian dynamics for most practical purposes and that the invariant manifolds 
constructed from them are "almost invariant" for the full Hamiltonian dynamics. 

We emphasise again that in this section we omitted superscripts (A) on the 
coordinates in order to keep the notation simple and that the integrals of the motion 
of the A**^ order normal form only assume the simple form in (j2.78|) if they are 
expressed in terms of the A**^ order normal form coordinates p.59p . 



Fi gure 1: The left figure shows contourhnes of the action angle variables / and (^i (hyperbolae and 
straight lines, respectively) in the saddle plane with coordinates {qi,pi) and {Qi,Pi) which are rotated 
versus each other by 45° degrees. The right figure shows contourlines of the action angle variables Jk and 
ipk, k = 2, . . . ,d, (circles and straight lines, respectively) in the centre planes with coordinates {qk,Pk) = 

{Qk,Pk)- 

3 Quantum Normal Form Theory 

In this section we develop a normal form theory for quantum mechanics that is 
algorithmically the same as the one presented for classical mechanics in the previous 
section. Sect. [2j However, the objects manipulated by the algorithm in the quantum 
mechanical case are different, and this is what we now describe. 

In quantum mechanics the role of a Hamilton function in classical mechanics is 
played by a self-adjoint operator, the Hamilton operator. While the Hamilton 
function in classical mechanics acts on a phase space, which was R^'^ in Sect. [21 a 
Hamilton operator acts on a Hilbert space, which will be L'^iM!^) in our case. 

The quantum mechanical analogue of a symplectic transformation in classical 
mechanics is a unitary transformation. The conjugation of a Hamilton operator H 
by a unitary operator U gives the new operator 

H' = U*HU, (3.1) 

where U* denotes the adjoint of U. The operator H' is again self-adjoint and has the 
same spectral properties as the original Hamilton operator H. We will use unitary 
transformations to simplify the Hamilton operator in the same way that we used 
symplectic transformations to simplify the classical Hamilton function. In the classi- 
cal setting the symplectic transformations were obtained as the time-one maps of a 
Hamiltonian flow, where the Hamiltonian, W, was referred to as the generating func- 
tion. In the quantum mechanical setting we will analogously consider a self-adjoint 
operator W which gives the unitary operator 

U = e-^^. (3.2) 

The operator W is called the generator of U. Analogous to the development of ()2.7p 
and the results that follow, we now consider the one parameter family of self-adjoint 
operators defined by 

#(e) :=e^^#e-^^, (3.3) 



where the parameter e is real. Note that H' = U*HU = H(e = 1), and H = H{e = 
0). If we differentiate (j3.3|) with respect to e we obtain the Heisenberg equation 

^H{e)='-\^,H{e)], (3.4) 

where [•,•] denotes the commutator which, for two operators A, is defined as 
[A^ B] = AB — BA. Therefore H' can be obtained from the solution of ()3.4p with 
initial condition H{e = 0) = H. Equation (j3.4p will play the same role for the 
development of the quantum normal form as Equation (12. 9p played for the classical 
normal form. This is consistent with the usual quantum-classical correspondence 
where the commutator ^[-j-] is related to the Poisson bracket {•,•}. In the next 
section we will make this correspondence more precise. 

One of the key properties of the classical normal form in the neighbourhood of a 
nonresonant saddle-centre-... -centre equilibrium point is that the Hamilton function 
in normal form is a function of the classical integrals, see (|2.83p . We will see in 
Sec. m that this feature will help us to understand the local classical dynamics and 
identify the phase space structures that control the dynamics near a non-resonant 
saddle-centre-... -centre equilibrium point. In the quantum mechanical case the clas- 
sical integrals will correspond to "elementary" operators with well known spectral 
properties. Analogously to symplectic transformations in the classical case, we will 
use unitary transformations in the quantum mechanical case to bring the Hamilton 
operator into a simpler form in which it will be a function of these elementary op- 
erators only. In the same manner as in the classical case, this simplification will 
be obtained "order by order". To give notions like 'order' and 'equilibrium point' a 
meaning for quantum operators and also to derive an explicit algorithm to achieve 
the desired simplification we will have to relate quantum operators to classical phase 
space functions and vice versa. This is the subject of the following section. Sec. 13.11 
The formalism developed in Sec. 13.11 is then used in Sec. 13.21 to transform Hamilton 
operators through conjugation by unitary operators. In Sec. 13.31 we will define when 
a Hamilton operator is in quantum normal form, and show how a given Hamilton op- 
erator can be transformed to quantum normal form to any desired order. In Sec. 13.41 
we study the nature of the quantum normal form for our case of interest, which is 
in a neighborhood of a non-resonant saddle-centre-- • • -centre equilibrium point of a 
corresponding classical Hamiltonian system. As a first explicit example, we show how 
the quantum normal form can be computed for one-dimensional potential barriers in 
Sec. [331 

3.1 The Classical- Quantum Correspondence 

The basis for our quantisation of the classical normal form described in Section 12.21 
is the Weyl quantisation and the associated Weyl calculus. Before we use the Weyl 
calculus to define the quantum normal form in Sec. 13.31 we want to give some back- 
ground on the general theory, which provides a systematic way of formulating the 
quantum-classical correspondence. General references for the material in this section 
that contain much more detail and background are |Fol89[ IDS991 IMar02] . 

3.1.1 Weyl quantisation 

A quantisation is a rule which associates operators on a Hilbert space with functions 
on a phase space. We will use here the Weyl quantisation, which is the one most 
commonly used. Let qk and pk, k = 1, . . . , d, he the components of the position and 



momentum vectors q and p, respectively. These are quantised in such a way that 
they act on a wavefunction ip{q) according to 



= qki^iq) , Pki^iq) = — • (3-5) 

The Weyl quantisation extends these prescriptions to general functions of {q,p) by 
requiring that, for ^q,^p G M!^, the quantisation of the exponential function 

es((«P-9>+(5<3.P» (3.6) 

is the phase space translation operator 

f^^^^^ = es(«f'9)+(€9.P)) . (3.7) 

Using Fourier inversion we can represent a function on phase space as 

A{q,p) = j:r^[ [ :4(e„e,)ei((«p.'?)+(««.p))de,dep, (3.8) 

where 

M<i,^p)= [ I ^(<?,p)e-i(«-'?>+«-^'>) dgdp (3.9) 

is the Fourier transform of A. The Weyl quantisation Op [A] of A is then defined by 
replacing the factor eft ^^^^"'^^^^^"^^^ in the integral (13. Sp by the operator T^g,^p, i.e. 

In order to manipulate these operators and understand their mathematical properties 
we will need the appropriate definitions and notation. We will say that A G Sfi(W^ x 
Mf^) if A depends smoothly on {h, q,p) and if for all a, /9 G M'^ and G IN there exists 
a constant Ca,/3,k such that 

(1 + + \p\)''\d^d^A{h,q,p)\ < . (3.11) 

The space Sfi{M!^ x M!^) is similar to the usual Schwartz space. The only difference is 
that we allow the functions to depend additionally on the parameter h in a smooth 
way. For A G Sfi(R!^ x R"^) the Fourier transform is again a Schwartz function and 
so the Weyl quantisation (j3.1Up gives a well defined bounded operator. But the 
quantisation can be extended to larger classes of functions. One such larger standard 
class of functions for which the Weyl quantisation is well behaved is S"^{M,'^ x M,'^) 
for some m G IR. Here A G S""(Il'^ x R*^) if A satisfies the estimates 

\d^d^A{h,q,p)\ < CaA^ + \q\ + \p\r for all a,/3 G M'^ . (3.12) 

If ^ G S"'(R'^ X R'^) then Op[A] : Sn{R'^) 5n(E'^) (see, e.g., [DH99]). Here 5n(E'^) 
is defined analogously to 5^(111'^ x R*^) in ()3.1ip . The function A is called the (Weyl) 
symbol of the operator Op[^]. If the symbol A also depends on the parameter h 
we will usually assume that, for small h, A has an asymptotic expansion in integer 
powers of h, 

A{h,q,p) ^ Ao{q,p) + hAi{q,p) + h'^A2{q,p) + . . . . (3.13) 

Here the leading order term Ao{q,p) is then called the principal symbol and it is 
interpreted as the classical phase space function corresponding to Op[A]. 



The quantisation ()3.10p can also be inverted. Let A be an operator, then 



A{h,q,p):=Tr{f*{q,p)A), (3.14) 

is the Weyl symbol of A, i.e. we have Op[^] = A, with Tr denoting the trace, and 
T* denoting the adjoint of T. 

The advantage of this representation of operators is that many properties of the 
operators are nicely reflected in their symbols. For later reference we collect two such 
relations: 

1. For the adjoint operator one has Op[A]* = Op[A*], where A* denotes the 
complex conjugate symbol of A. Hence, a real valued symbol gives a symmetric 
operator. 

2. If A € S'^(]R°' X ]R°'), i.e., the symbol and all its derivatives are bounded, then the 
corresponding operator is bounded as an operator on L^(]R'^). This is known 
as the Calderon-Vaillancourt theorem [DS99j . This implies in particular that a 
real valued symbol A G 5'^(1R'^ x M!^) gives a self-adjoint operator Op[^]. 

For example, the symbol J = {p^ + q^) /2 on IR x R is in ^^(IR x IR). Its principal 
symbol is (p^ + q^^jl and the Weyl quantisation gives 

Similarly, the symbol / = pt; is in S'^(l[l x H) with principal symbol pq and is 
quantised as 

n( d 1^ 



These are the quantisations of the classical integrals obtained in Section 12.3.21 and 
they will form the building blocks of the quantum normal form associated with a 
saddle-centre-- • • -centre equilibrium point in Section [331 



3.1.2 The Moyal bracket 

The main idea behind the introduction of symbols of operators is that one can use 
the symbols to study properties of the operators, as we already indicated in the 
last subsection. Since the symbols are functions they are in general much easier 
to study than operators. One can probably say that the single most useful fact 
about pseudodifferential operators, i.e., operators whose symbols satisfy estimates 
like ()3.12p . is that they form an algebra, i.e., the product of two such operators is 
again of this type, and that one can compute the symbol of a product from the 
symbols of the operators which are multiplied. 

The quantum normal form algorithm we will develop will rely essentially on this 
product formula for symbols. Given two functions j4, i3, one can find a function A*B 
such that Op[A] Op[i?] = Op[A*i?], see [DS99j . This so called star product of A and 
B is given by 

A * B{q,p) = A{q,p) exp (^i^[{dq, dp) - {dg, dp)]^ B{q,p) , (3.17) 

where the arrows indicate whether the partial differentiation acts to the left (on A) or 
to the right (on B). For the precise meaning of the expression on the right hand side 
of this equation we refer the reader to [Fol89t[DS991 IMar02j . However, by expanding 



the exponential we obtain the more exphcit asymptotic expansion in powers of h that 
will suffice for our purposes 



A * B{q,p) ~ ^ - i- A{q,p)[{dg, dp) - {dg, dp)fB{q,p) 

k=o "-^^^ (3.18) 

= A{q,p)B{q,p) + ^{A, B}{q,p) +■■■ , 

where {•,•} again denotes the Poisson bracket defined in (|2.8|) . In particular, if 
A G 5™(E'^ X R'^') and B G S""' {R'^ x E'^) then A*B e 8""+""' (R'^ x R'^) f [Ebl89l 
IDS99[ IMar02j ). It is worth mentioning that even if A and B are independent of h, 
the product A* B will in general depend on h with the principal symbol being given 
by AB, i.e., the usual product of the functions A and B. 

From the Heisenberg equation p.4p we see that the commutator plays an impor- 
tant role when one wants to conjugate an operator with a one-parameter family of 
unitary operators. Applying the product formula (j3.17p to the expression for the 
commutator of Op[^] and Op[i?], 

Op[^] Op[B] - Op[B] Op[A] = Op[A *B]- Op[B * A] = Op[A * B - B * A] , (3.19) 

we obtain the formula for the symbol of a commutator 

{A*B-B*A)iq,p) = ^ {A, S}^ {q,p), (3.20) 

where {■,-}m is the Moyal bracket which is defined as 

{A,B]M{q,p) ='^A{q,p)s\nl^[{^p,1)q) - {dp,^q)\^B{q,p) . (3.21) 

For the precise interpretation of the right hand side of this equation we again refer 
the reader to ^Fol89t [DS991 [Mar02j . However, as above, by expanding the sine we can 
obtain an explicit asymptotic expansion for small h that will suffice for our purposes, 

{A,B}M{q.p) - [2) (^^Tiyi'^^'^'^^f^^^' ~ d,)\^^'^'^B{q,p) . 

(3.22) 

Note that in the case where one of the functions A,B is a polynomial the sum 
terminates at some finite k and gives the exact expression for the Moyal product. In 
what follows, all our explicit calculations will use from the Weyl quantisation only 
the asymptotic formula (j3.22p for the Moyal product. Since we will only work with 
finite Taylor series, the asymptotic expansion will always terminate and give the 
exact result. 

From (j3.22p we see that 

{A,B]M{q,p) = {A,B]{q,p) + Oih") , (3.23) 

i.e., in leading order the Moyal bracket is equal to the Poisson bracket, and moreover, 
if at least one of the functions A,B \s a, second order polynomial then 

{A,B}M{q,v) = {A,B}{q,p) . (3.24) 



3.1.3 Localising in phase space 



One important application of the product formula ()3.17l) is that it allows to localise 
operators in phase space, a technique often called micro-localisation, which in fact 
gave the whole field of microfocal analysis its name. Let p G 5°(R'^ X R"^) be a cutoff 
function, i.e., there is a set U C R*^ x R'^ such that 



and p has support in a small neighbourhood of U. Then we will call Op[/3] a cutoff 
operator (associated with U), and we can use it to split any operator Op[-fr] into two 
parts 



Op[H] = Op[p] Op[H] + (1 - Op[p]) Op[H] = Op[Fioc] + Op[H,em\ (3-26) 



where -ffioc = p * H and //rem = H — p * H . By the product formula (j3.17p the 
symbol Hioc is concentrated near the support of p and //rem is concentrated on the 
complement of the support of p. In this sense the usual procedure to localise the 
study of functions and dynamical systems by multiplication with cutoff functions 
can be quantised. In particular we have -fTioc = pH + 0{h), so the leading order is 
actually the classical localisation. If Op[p] is a cutoff operator associated with some 
phase space region U we will call Op[/fioc] = Op[p] Op[H] the localisation oi H to U. 

The localisation appears to be a very natural object to consider with regard to 
the application we are interested in, namely the study the dynamics of a chemical 
reaction which is described by a Hamilton operator Op[-ff] whose principal symbol 
has a saddle-centre-- • • -centre equilibrium point. The neighbourhood of the equilib- 
rium point is the most important region for the chemical reactions. This is where 
the reactants combine to form the activated complex which then decays into the 
products. So it is natural to use the above procedure to localise the Hamiltonian 
to a neighbourhood of the equilibrium point in phase space. In fact, we will derive 
the quantisation of the classical normal form procedure for a Hamiltonian which is 
localised. 

The localisation has another advantage which is of a more technical nature. The 
Hamilton operators we will encounter have symbols with polynomial growth in p and q 
for large p and q, and this leads to some technical complications concerning questions 
like self-adjointness and unitarity. If we localise our Hamiltonians by multiplication 
with a cut-off operator we end up working with operators with bounded symbols only, 
for which self-adjointness is easy to show. This will make many proofs technically 
much easier. 

3.2 Transformation of Operators through Conjugation 
with Unitary Operators Using the Weyl Calculus 

We will now apply the Weyl calculus to the problem outlined in the beginning of this 
section. For an operator A = Op[A\ with symbol A we consider its conjugation by 

a unitary operator U = e^^, where W = Op[W] has symbol W. Our aim is to find 
the symbol A' such that 



P\u = l 



(3.25) 



Op[A'] = et«P[^10p[A]e-iOp[^l . 



(3.27) 



If we introduce the one parameter family of operators 



A{e) = Op[A{e)] = e^^Pl^] Op[A]e-r^'^Pl^^ 



(3.28) 



then A' = A{€ = 1) and the Heisenberg equation ()3.4p can be written in terms of the 
Moyal bracket as an equation for the symbol ^(e), 

^A{e)={W,A{e)}M ■ (3-29) 
de 

In order to obtain A' we thus have to solve (I3.29P with initial condition A{<S) = 
A. Note the similarity between (j3.29|) and (j2.9|) in Section 12.11 which expresses the 
correspondence between the Heisenberg equation ()3.4p and the classical equation 
(j2.9p in the framework of the Weyl calculus. 

We will now discuss methods of how to solve equation (j3.29p for certain choices 
of W . Recall that if is a polynomial of order less than or equal to two, then the 
Moyal bracket reduces to the Poisson bracket (see (|3.24p ) and hence Equation (j3.29p 
reduces to (12. 9p . and we recalled earlier in our development of the classical normal 
form theory that polynomials of order less than or equal to two generate affine linear 
symplectic transformations (see Sec. l2.2l and reference [Fol89] ) . The following Lemma 
tells us that the symbols of operators transform in the same way as classical phase 
space functions under such transformations. 

Lemma 6 (Exact Egorov). Assume W{q,p) is a polynomial of order less than or 
equal to 2 with real valued coefficients, and let be the time one map of the Hamil- 
tonian flow generated by W (see (|2.5p ). Then 

is unitary, and for every A G S"^(R!^ x W^), we have 

Op[^l Op[A]e-i °P[^1 = Op[A'] (3.31) 
with A' E 5"'(]R^ X R'^) given by 

A' = Ao'^:^ . (3.32) 

Proof. For the full proof we refer the reader to the appendix to chapter 7 in [DS99] . 
The main ideas are as follows. It is well known that Op[IF] is essentially self-adjoint 
(see e.g. [DS99]), and therefore U is unitary. In order to find H' we have to solve 
(j3.29p . Since is a polynomial of order two or less than two Equation (j3.29p reduces 
to ()2.9p . From Equation ()2.5p we see that A{e) = Ao and at e = 1 we obtain 
(j3.32p . Now if is a polynomial of order less than or equal to two, then is an 
affine linear transformation. Hence if ^ G 5™(]R'^ x R'^), then A' G S'™(]R'^ x R'^). □ 

This result is called 'exact Egorov' because there is a more general theorem due 
to Egorov |Ego69| which states that, for a large class of W, a similar result holds 
asymptotically for ^ ^ 0. However, only for polynomials of degree equal to or less 
than two, the higher order terms in h vanish. 

For later reference we consider the following example. For {q,p) £ R^, let 



W{q,p) = -^^{p^+q^), (3.33) 

which is the Hamilton function of an harmonic oscillator. The factor — 7r/4 is intro- 
duced for convenience. The function W generates the vector field 

Xw{q,p) = {dpW{q,p),-dgW{q,p)) = ^{-p,q). (3.34) 



The corresponding flow is given by 

IT IT IT IT 

(9(e),p(e)) = '^wi.Q^P) = (cos(e-)g-sin(e-)p,sin(e-)g + cos(e-)p) . (3.35) 

The harmonic oscillator thus generates rotations in the (g,p)-plane. In particular, 
the time one map of the flow generated by W gives the map from the coordinates 
{q,p) to the new coordinates 

{Q,P) = ^Uq,p) = ^{Q-P,Q + p). (3-36) 

which we already considered in Sec. 12. 3i Transforming I{q,p) = pq under this flow 
we get 

I\Q, P)=Io $^i(Q, P) = i(p2 _ q2^ ^ (3.37) 



which gives the operator 



We will refer to Op[/'] as the Q representation of Op[/], and for later reference we 
denote the unitary transformation which classically generates the 45° rotation (j3.36p 
as 

= e-^Op[W^], (3.39) 

where W is given by (j3.17p . 

Having discussed this particular example of an application of Lemma [6] (exact 
Egorov) we now turn to the case of higher order polynomials in W. To this end we 
will develop a power series approach analogous to what we described in Sec. 12.11 This 
will provide higher order approximations in h as well as give an explicit expression for 
the symbol of the transformed in ()3.27p in the case where or j4 are polynomials. 

We begin by simplifying the notation and define the Moyal-adjoint action. For 
two smooth functions W and A, we define analogously to the adjoint action in p.lOp 
the Moyal-adjoint action as 

Madw A :={]¥, A}m ■ (3.40) 
Using the Moyal adjoint Equation (j3.29p becomes 

dA{e) 



de 



MadwA{e). (3.41) 



We compute higher order derivatives of A(e) with respect to e in a manner analogous 
to (12. lip and ()2.12p . We successively differentiate ()3.29p and apply the notation 
([OO]) to obtain 

^A{e) = [Ms.dwrA{e) . (3.42) 
Hence, the (formal) Taylor series in e around e = is given by 



A{e) = Y — \MadwrA , (3.43) 
and setting e = 1 we obtain the formal sum 

A' = ^ — [Madvi/]"^ . (3.44) 

n=0 ^' 



This expression is completely analogous to ()2.13p . and as we will see in more detail, 
can be used in a similar fashion to compute the symbol A' up to any desired order in 
h and {q,p). In particular, analogous to equation ()2.13p . it gives the Taylor expansion 
with respect to e, evaluated at e = 1, for the symbol A' of the operator obtained after 
conjugation of the operator defined by the symbol A by the unitary transformation 
generated by W. This formula forms the basis of the quantum normal form method 
where the idea is to "simplify" (or "normalise") the symbol whose quantisation will 
then correspond to the normal form of the Hamilton operator. As in the classical 
case, the computation of the Taylor expansion is carried out "order by order" using 
power series expansions of the symbol in {q,p) and h. The series is expanded about 
an equilibrium point of the principal symbol, and therefore the quantum normal 
form will be valid in a neighbourhood of this point. Hence, as in the classical case, 
the quantum normal form is a "local object" whose operator nature requires more 
technical details for a rigorous characterization of its properties (cf. Section 13.1.31 
and Definition U]) , and we will describe these in more detail in the following. 

Therefore similar to the mathematical formalism required for computing the clas- 
sical normal form, normalising the symbol of the operator that will correspond to 
the quantum normal form will require us to manipulate monomials which in addition 
to {q,p) now also have factors of h. In order to describe this we adopt a notation 
introduced by Crehan [Cre90] and define the spaces 

= span U^q^p^ := W J] pf*^ : \a\ + + 2j = s\ . (3.45) 
fe=i ^ 

These spaces Wqjjj are closely related to the spaces W^fj spanned by the polynomials 
()2.14p in the classical case. In fact we have 

[s/2] 
fe=0 

where [s/2] denotes the integer part of s/2. 

Below we want to use functions W £ Wqjn in order to construct unitary operators 

of the form e~ ft ^pI^I. However, the quantisation of a function W £ Wqj^ will give an 
unbounded operator and this makes the discussion of self-adjointness of Op[Ty], and 
hence the unitarity of e~ ft ^P'^^l, more complicated. But since we are interested in the 
local quantum dynamics generated by a Hamilton operator in the neighbourhood of 
an equilibrium point of its principal symbol it will be sufficient to have a local version 
of the spaces Wq^. We thus apply the localisation procedure from Section [3.1.31 We 
say that W G Wqj^. j^^, ifWe Sfi(^'^ x K.'^) and there is an open neighbourhood U of 
zo = eR'^ xR'^ such that 

W\u e . (3.47) 
The quantisation of elements of Wq^^. j^^. will then give bounded operators. Therefore, 

ifWG VVqm.ioc is ^^^^ valued then Op[l^] will be self-adjoint and thus U = e~Ti'~'P^^^ 
will be unitary. 

We will frequently use Taylor expansions and want to modify them in such a way 
that the terms in the expansion are in Wq^^. j^^. In order to make our discussion of 
this property precise we will need the following definition. 

Definition 4. We will say a function On G ShiR!^ x W^) is a remainder of order 
N (around {q,p) = (0,0) j if there is an open neighbourhood U of {q,p) = (0,0) and 



c > such that 

\ON{e'^h,eq,ep)\ < ce^ (3.48) 

for h < I, {q,p) G U and e < 1. 

We then can formulate 

Lemma 7. Let A G Sfii]^'^ x ^'^), then there exist As G VVqm-ioc ^'^'^^ ^^^'^ f^''" '^^2/ 
N & there is a remainder On G •S^i^'^ x R*^) of order N such that 

N-l 

A=^As + On- (3.49) 

s=0 

Proof. Let us take the ordinary Taylor expansion of A{h,q,p) around {h,q,p) = 
(0, 0, 0) and order the terms according to the definition of order in (|3.45|) . This gives 
us an expansion A = Xlfio^ As + Rn with 

^« = E 1^^^^^"^^^^°' ^o,Po)qy^h^ e W^n, (3.50) 

and R]\!{e'^h,eq,ep) = 0{e^). Now choose a function p G 5^(11'^ x R'^) with p|[/ = 1 
for some open neighbourhood U of 0, and set Ag := pAs- Then it follows directly 
that As G yV^^.ioc and O^v := A - ^l^'o^ S 5a(]R'^ x R'^) is a remainder of order 
A^. ' □ 

The main reason for defining the order s according to (j3.45p . i.e., the reason for 
double counting the powers of h, is that it behaves nicely with respect to the Moyal 
product. This is reflected in the following lemmata. The first one is the analogue of 
Lemma [T] in the classical case. 

Lemma 8. Let W G >V^;„^i„,, A G W^^^i^,, s,s' > 1, then 

{Ty,A}MGW^+f;-2, (3.51) 

and for n > 0, 

[MadH.]"^G>vi;;?+% (3.52) 
if n{s' — 2) + s > and [Mad^y = otherwise. 
Proof. We can write the Moyal bracket (j3.22p as 

{W,A}m = E )''(^^^^''^'^(^'^)(«'^) (3.53) 
with the bi-differential operators 

D^''^+'\W,A){q,p) := Wiq,p)[{dp,d,) - {dp,d,)f'^'^A{q,p) , (3.54) 

Now the bi-differential operator jg q£ Qj-^^^gj- 2k + 1 in the arguments involving 

A and W individually, and therefore 

D^''^'^ : W^^^,„, X W(^.,^^ ^ K^'o:''^'-''"^'' ■ (3.55) 

On the other hand, multiplication by h^^ maps y^^-^'io^^^^^ {2fc+i) 
^qm*^?oc^^^^* (2fc+i)+4fc _ yys+s^-2^ therefore every term in the series (13.530 is 



in VVq^!*[~^^. But the order of W and A as polynomials in {q,p) near {q,p) = (0,0) 
is at most s and s', respectively, and therefore the terms in the series (|3.53|) vanish 
near {q,p) = (0,0) for 2k + 1 > min(s,s'). Hence 

W^lAf ew^ifc'- (3-56) 
The second result follows then by induction. □ 

We can now turn our attention to the computation of the symbol of a conjugated 
operator when the generator of the unitary operator has order larger than 2. The 
computation will proceed in two steps, in the first lemma we show that conjugation 
respects the class of symbols we are working with. 

Lemma 9. Let W G y^qm ioc '^'^'^ ^ ^ ■Sni^'^ x R"'), then there exists an A' £ 
SniM,"^ X R'^) such that Op[A'] = eT^^^^^^ Op[A]e-i°P['^l. 

The techniques for proving this lemma are different from the ones we use in the 
rest of the paper. In order not to interrupt the flow of the paper, we therefore present 
the proof in Appendix Rl 

By Lemma [9] we know that the symbol of e^ ^p^^^' Op[A]e~ « ^p''^' is a function 
in S!i(R!^ X R*^). With the help of Lemma [8] we can reorder the terms in the formal 
expansion (|3.44|) to turn it into a well defined Taylor expansion in the sense of Lemma 
[71 This is the content of the following Lemma which can be considered to be the 
analogue of Lemma E] in the classical case. 

Lemma 10. Let W G y^qai-loc ^' — 3; ^^^^ ^ ^ '5a(R'^ x R'^) with Taylor expansion 
A = 127=qAs, As e Then the symbol A' o/e^Op'^l Op[^]e-i OpI^l has the 

Taylor expansion 

oo 

^' = E ^'^ (3.57) 

with 

= Yl ;^[Madw^]'^^_„(y„2) e Wqm;loc , (3-58) 
n=0 

i.e., for every N G H there exists a remainder Oat G 5ft(R°' x R'^) of order N such 
that 

N-l 

A' =Y^ ^'s + On ■ (3.59) 

Proof. By Lemma [9] we know that A' € 5?i(R'^ x R*^), and we have to compute its 
Taylor series. With (j3.42p we can use the Taylor expansion of A'{e) to write 



A' =J2^ [ M&dw ] + O'j^ (3.60) 



n=0 

with 



O'n = (^^^ - e)^-'[MadwfA'{e) de , (3.61) 

being just the standard remainder formula for Taylor expansions. Since ^'(e) G 
5ft(R'^ X R'^) = we have by Lemma [H] that G Sn(R'^ x R'^) is a remainder 

of order iV. If we next insert the Taylor expansion for A we get 



N-lN-l ^ 



A' =}^}_^-[MadwrAi + ON, (3.62) 

;=o n=0 



where Oat G Sfi{W^ x M!^) denotes the cohection of all the remainder terms of order 
N. Using Lemma [5] we can collect all the terms of order k in the sum which gives 
(I3.58p . To this end one can proceed completely analogously to the proof of Lemma [2] 
and we therefore omit the details. 

□ 

3.3 Definition and Computation of the Quantum Nor- 
mal Form 

We will now define when a Hamilton operator is in quantum normal form. Similar 
to the case of the classical normal form, in general a Hamilton operator is not in 
quantum normal form. However, as we will show, the formalism based on the Weyl 
calculus developed in the previous two sections can be used to construct an explicit 
algorithm which will allow us to transform a Hamilton operator to normal form to 
any desired order of its symbol. The algorithm will consist of two parts. The first 
part operates on the level of the symbols of operators, and this part of the algorithm 
will be very similar to the normalisation algorithm in the classical case. In the second 
part the symbols are quantised, i.e., the operators corresponding to the symbols will 
be determined. 

The starting point is a Hamilton operator Op[i?] which is the Weyl quantisation 
of a symbol H(h,q,p). Assume that the Hamiltonian dynamical system defined by 
the principal symbol has an equilibrium point at zq = {qo,Po), i-e., the gradient of the 
principal symbol vanishes at zq. Let H2{z) G Wqj„ denote the second order term of 
the Taylor expansion of the symbol H about zq and Op[/72] its Weyl quantisiation. 
We now make the 

Definition 5 (Quantum Normal Form). We say that Op[if] is in quantum normal 
form with respect to the equilibrium point zq of its principal symbol if 

[Op[i72],Op[//]] =0, (3.63) 

or equivalently in terms of the symbol, 

adH^ H = {H2, H} = . (3.64) 

The equivalence of the two equations in Definition [5] derives from the fact that 
the Moyal bracket reduces to the Poisson bracket if one of its argument is quadratic. 
Moreover, we remark that H2 and the quadratic part of the principal symbol differ at 
most by a term that consists of h with a constant prefactor. Since the Poisson bracket 
vanishes if one of its two arguments is a constant it does not make a difference in 
Definition [5] if H2 in (j3.64p would be replaced by the second order term of the Taylor 
expansion of the principal symbol. 

Like in the case of a Hamilton function being in classical normal form the property 
of a Hamilton operator to be in quantum normal form has strong implications which 
in the quantum case lead to a considerable simplification of the study of the spectral 
properties of the operator. To this end recall that two commuting operators have 
a joint set of eigenfunctions. Hence, if an operator is in quantum normal form the 
study of its spectral properties will be simplified considerably, since the spectrum 
and eigenfunctions of an operator Op[/72] with a symbol of order 2 are well known. 

Similar to the classical case a Hamilton operator is in general not in quantum 
normal form. However, we will now show how the formalism developed in the previous 
two sections can be used to transform a Hamilton operator to quantum normal form 



to any desired order of its symbol. Similar to the classical case we will truncate 
the symbol at a certain order and show that the corresponding Hamilton operator 
will lead to a very good approximation of many interesting spectral properties of the 
original Hamilton operator. 

We develop the following procedure. Let H = H^'^'> denote the symbol of our orig- 
inal Hamilton operator. We will construct a consecutive sequence of transformations 
of the symbol according to 

H =: ^ ^ ^ ^(3) ^ . . . ^ H^N) ^3_g5) 

by requiring the symbol H^'^\ for n > 1, to derive from the symbol by 
conjugating Op[H^'"'~^^ with a unitary transformation according to 

Op[F(")] = e^ °P[^"1 Op[F("-^)]e-^ 0P['^"1 , (3.66) 

where the symbol Wn of the generator of the unitary transformation is in JA^. j^^,. 
Like in the series of symplectic transformations in the classical case in (I2.27p . N 
in (j3.65p is again a sufficiently large integer at which we will truncate the quantum 
normal form computation. The algorithm for normalising the symbol will be identical 
to the classical case. The key difference is that the Poisson bracket of the classical 
case is replaced by the Moyal bracket in the quantum case. With this replacement, 
the mathematical manipulations leading to normalisation of the symbol are virtually 
identical. 

Towards this end, using ()3.44p we see that, analogously to ()2.35p in the classical 
case, we have 

oo ^ 

^^"^ = E ^ [Madty„ . (3.67) 

fc=0 

Like in the classical case the first two steps, n = 1,2, in (I3.65|) differ somewhat 
in nature from the steps for n > 3. The first step serves to shift the equilibrium 
point to the origin and the second step serves to simplify the quadratic part of the 
symbol. It follows from Lemma (exact Egorov) that we achieve these affine linear 
transformations by choosing the symbols Wi and W2 identical to the generators of 
the corresponding symplectic transformations in the classical case. We thus have 

H^^\h,z) = H^^\h,z + zo) (3.68) 

and 

H'^'^\h,z) = H^^\h,M-^z), (3.69) 

where M is a suitable symplectic 2d x 2d matrix which achieves the simplification of 
the quadratic part of the symbol analogously to the classical case. It is important to 
note that we do not explicitly need the generators Wi and W2 which, as mentioned 
in Sec. 12. 2^ might be difficult to compute. 

Before we proceed with the normalisation of the higher order terms, n > 3, we 
will assume that we localise around the equilibrium point which is now at the origin, 
see Section [3. 1.3t i.e., by multiplying //^^^ by a suitable cutoff function concentrated 
about the origin we can assume H^"^^ G ShiJ^'^ x m*^) and the terms i^i^^ of the Taylor 
expansion of H^^^l to be in W^j^. j^^. 

For the higher order terms, n > 3, we find by ()3.58p in Lemma [TOl that the terms 
Hs""^ can be computed from the terms of the power series of according to 

H^""^ = E ^[Mad^J^i^t^d,) . (3.70) 

k=0 



The normalisation procedure for the terms of order n > 3 of the symbol has very 
similar properties as the corresponding procedure in the classical case. In particular a 
transformation at a given order does not affect lower order terms. This is made more 
precise in the following lemmata that are the analogues of Lemma [3] and Lemma H] 
for the classical case from Section 12.21 

Lemma 11. = iff \ n > 3. 

Proof. The proof is completely analogous to the proof of Lemma O and is therefore 
omitted. □ 

Like in the classical case, Lemma [H] motivates the adoption of the following 
notation for the operator 

P:=ad^(2) ={iif ,•}. (3.71) 
Lemma 12. For n> 3 and < s < n, Hi"^ = H^" ^\ 

Proof. The proof is completely analogous to the proof of Lemma S] and is therefore 
omitted. □ 

Like in the classical case the n*^ order term in indicates how to choose Wn 
for n > 3. 

Lemma 13 (Quantum Homological Equation). For s = n > 3, 

Proof. The proof is completely analogous to the proof of Lemma [5] and is therefore 
omitted. 

□ 

The homological equation ()3.72p is solved in exactly the same way as the homo- 
logical equation in the classical normal form computation described in Sec. 12.3. 1[ 
The only difference is that we now deal with a symbol that in contrast to the clas- 
sical Hamilton function in general depends on h. But due to the splitting = 

©Jf^Q /I'^W^j"^'^, see (|3.46p . the results on the solution of the classical homological 
equation can be transferred directly. In particular the notion of solvability introduced 
in Definition [2] carries over verbatim. 

We note that so far we have only shown how to transform the Hamilton operator 
to quantum normal form on the level of its symbol. We have not yet discussed the 
implications for the corresponding transformed operator. As we will see, similar to 
the question of how to explicitly solve the homological equation, the nature of the 
transformed Hamilton operator depends on the type of the equilibrium point of the 
principal symbol. In the next section, Sec. 13.41 we will discuss this in detail for the 
case of a non-resonant saddle-centre-- • • -centre equilibrium point. 

We summarise our findings in the following 

Theorem 2. Assume the principal symbol of Op[H] has an equilibrium point at 
zq S M!^ X M!^, and that the homological equation is solvable in the sense of Def. [3 
Then for every € IN there is a unitary transformation Un such that 

U*^ Op[H]Un = OpIh'^qJp] + Op[Oiv+i] (3.73) 



where OplHqJp] is in quantum normal form (with respect to 0) and On+i is of order 
N + 1. 



Proof. As we have seen in this section the conjugations of a Hamilton operator by 
unitary transformations to transform it to quantum normal form can be carried 
out on the level of the symbols of the operators involved. This makes the proof of 
Theorem [2] very similar to the proof of Theorem [T] in the classical case. In fact, 
the proof of Theorem [T] carries over verbatim when one replaces the Poisson bracket 
by the Moyal bracket. Then Lemma [3] is replaced by Lemma [TT] and Lemma H] by 
Lemma [121 

Using the scheme (|3.65|) with (|3.66|) gives then the unitary transformation Un in 
(frajl as 

^^ = e-tOp[m]e-tOp[T4^2]g-iOp[Ty3]...e-iOp[Tyiv] . (3.74) 

The first two generators, Wi and W2 are chosen exactly as in the classical normal 
form algorithm, see (j3.68p and the following paragraph, by Lemma [H] (exact Egorov) 
this induces the same transformation of the symbols as in the classical case. The other 
generators Wn, n > 3, are then chosen recursively as solutions of the homological 
equation, see Lemma[l3l where after each step we have to determine H^^^ up to order 
N from (137701) . 

□ 

Similar to the classical case the definition of the quantum normal form in Defin- 
tion [5] is of little value for practial purposes since we cannot expect the quantum 
normal computation to converge if we carry it out for N ^ 00 as required by Defin- 
tion [5l For applications it is more useful to consider the truncated quantum normal 
form. 

Definition 6 {N^'^ Order Quantum Normal Form). Consider a Hamilton operator 
Op[H] whose principal symbol has an equilibrium point at zq G x R"^ which, 
for G IN, we normalise according to Theorem [21 Then we refer to the operator 
Op[Hq2p] in Equation ()3.73p as the N^^ order quantum normal form (QNF) of 

We have seen that the procedure to construct the quantum normal form is very 
similar to the procedure to compute the classical normal form. In particular the 
homological equations (|2.46|) and (|3.72p which determine the choice of the succes- 
sive transformations p. 270 and (13.650 . respectively, look identical since the Poisson 
bracket reduces to the Moyal bracket if one of its argument is a polynomial of or- 
der less than or equal to 2. However, it is important to point out that this does 
not mean that the Moyal bracket completely disappears from the procedure in the 
quantum case. In fact, while the normalisation transformation at a given order does 
not modify lower order terms, it does modify all higher order terms, and the Moyal 
bracket plays an important role in this, see (I3.70p . Consequently the terms in the 
Taylor expansions of i/^") and the generators Wn will in general depend on h. 

Since the Moyal bracket tends to the Poisson bracket in the limit /i ^ we 
expect that the symbol of the quantum normal form should tend to the classical 
normal form, too. This is indeed the case. 

Proposition 1. The principal symbol of the N^^ order quantum normal Op[//q^^] 
is the classical normal form of order N , i.e., 



H'^^j)p{h,q,p) = H^^j^p{q,p)+0{h) . (3.75) 



Proof. This follows from an inspection of the construction of the classical and quan- 
tum normal forms. The first two steps are identical by Lemma [6] (exact Egorov). The 



homological equation determining the choices of the Wn is as well identical. What 
is different however is the transformation of the higher order terms, k > n. Here we 
have Equation ()2.36p in the classical case and Equation ()3.70p in the quantum case, 
and these equations differ by the use of the adjoint versus the Moyal adjoint. But 
since Madvi/ A = adw A + 0{h) and therefore 

Mad'^ A = ad'^ A + 0{h) (3.76) 

the differences in the higher order terms between the classical and the quantum trans- 
formation schemes are always of order h. This implies that the difference between 
the symbol of the quantum normal form and the classical normal form are of order 
h. □ 



3.4 Nature and Computation of the Quantum Normal 
Form in a Neighbourhood of an EquiUbrium Point of the 
Principal Symbol of Saddle- Centre-- • -Centre Type 

We now describe how the quantum normal form of a Hamilton operator can be 
computed in the case where the principal symbol has an equilibrium point of saddle- 
centre-- • • -centre type, i.e., the matrix associated with the linearisation of the Hamilo- 
nian vector field generated by the principal symbol has two real eigenvalues, ±A, and 
d — 1 complex conjugate pairs of imaginary eigenvalues iiwfc, k = 2, . . .d. We will 
assume that the uJk, k = 2, . . . d, are nonresonant in the sense that they are linearly in- 
dependent over the integers, i.e., k2UJ2 + - .. + kdUJd / for all (fcs, . . . kd) € ^'^-^-{0}. 

As mentioned in the previous section it follows from Lemma [6] (exact Egorov) 
that we can use the same affine linear symplectic transformations that we used in the 
classical case in Sec. [2] to shift the equilibrium point to the orgin of the coordinate 
system and to simplify the second order term of the symbol. We thus have 

oo 

i/(2)=i?o + Ff , (3.77) 

where 

d 

Hi'\h,q,p)=Xqm + Y,'^{pl+ql) + ch, (3.78) 

k=2 

where c is some real constant. 

We note that in terms of the coordinates {Q, P) we defined in Sec. 12.31 H2 is 
given by 

i/f {h,q,p) = \ (Pf - Q?) + E Y {Pi + Ql) + ' (3-79) 

k=2 

which is the analogue of Equation ()2.63p in the classical case. 
3.4.1 Solution of the homological equation 

We will solve the homological equation in the spaces Wqj^. The solution will then 
be localised by multiplication with a cutoff function afterwards to obtain elements in 
^qm ioc- This will ensure that the quantizations of these symbols are bounded and 
generate unitary operators. 



In order to solve the homological equation in Lemma [13] we perform the hnear 
symplectic complex change of coordinates {q,p) i-^ ix,^,) given by xi = qi, = pi 
and 

Xk := -^{qk - m) , ■= ^{Pk - iqk) , k = 2,...,d. (3.80) 

In terms of these coordinates the operator D defined in ()3.7ip assumes the simple 
form 

d 

V = A(^i% - xia^J +^1^^(49^^ - x^a^rj . (3.81) 

k=2 



In terms of these coordinates the spaces defined in (13.450 are given by 

= span [Wx^'i^ := xf^l" : \a\ + + 2j = n\ , (3.82) 

fc=i ^ 

and the operator T? acts on an element Wx°^^^ G according to 

^ n = (^(/^i - «i) + E - "'^)) n • (3-83) 

k=l ^ fe=2 ^ fc=l 

This means that the map T> can again be diagonalised and similar to the classical 
case we have that Wqj^ is given by the direct sum of the kernel of T) acting on Wqjj^, 
Ker I'l-^n and the image of V acting on W" , Im ^^|-^„ , i.e., 

W^^ = KerP| elmPj . (3.84) 

Now we can express H^^ as 

where H^.-^J E Ker and H^-jj^ E Im I^Ly^ . We can therefore choose Wn 

such that 

-^Wn = i^S"'^ , (3.86) 

and therefore 



Ht^=H^:Kel (3.87) 



Similar to the classical case the choice of Wn is not unique since one can always add 
terms from the kernel of . However, we will require Wn € Im ^'l-^n , i.e., we 

will invert V on its image Im ^^l-wyn . 

''^ ^qm 

Using our assumption that the frequencies U2, . . . ,LOd are nonresonant, i.e., lin- 
early independent over we see from (I3.83P that a monomial h^x°'^^ is mapped 
to zero if and only if Uk = (3k, k = 1, . . . ,d. In particular Ker 'Dly^^ = {0} if s 

is odd. This implies that unitary transformations can be constructed such that all 
odd order terms in the symbol of the conjugated Hamilton operator are eliminated. 
Moreover, for s even, the terms that cannot be eliminated are those which are sums 

(s) (s) 

of monomials for which xj^ and have equal integer exponents for all k = 1, . . . ,d. 

Concretely, we can compute Wn from (|3.86|) as follows. We assume that H^n-ira 
is the linear combination of L monomials of order n. 



<7:^=E/^/^^'n-r€"s (3.88) 

1=1 k=l 

with 2ji + Ylk=i + /3fc;i = ^ for all Z = 1, . . . , -L, and for all Z = 1, . . . , L, there 
is at least one = 1, . . . , d for which au-i ^ (ik;i (i-e., the vectors {ai-u • • • > oid-i) and 
. . . , /3rf;;) are different for all I = 1, . . . , L). Upon inspecting (j3.83p . and using 
(j3.86p . we see that a suitable generating function is given by 

^n=j: — — . — -n- n • (3-89) 

1=1 Kp1;1 - "1;/) + Lfc=2 i^fc(Pfc;« - Otk-i) 

As mentioned above this solution of the homological equation is unique if we 
require Wn to be in Im 'Dly.^ . 

3.4.2 Structure of the Hamilton operator in A^*"^ order quantum 
normal form 

In the previous section we have seen how to obtain the quantum normal form to order 
N in the case where the equilibrium point is of saddle-centre-- • • -centre type. So far 
these computations were carried out on the level of the symbols of the Hamilton 
operators. We now discuss the implications for the structure of the corresponding 
Hamilton operator in quantum normal form itself. 

In the classical case in Sections 12.3.11 and 12.3.21 we have shown that in each mono- 
mial of the Hamilton function in A^**^ order classical normal form the coordinate pairs 
{xk,^k) (or equivalently {qk,Pk)), k = 1,... ,d, occur with equal integer exponents 
and that this implies that the Hamilton function in A^*'^ order classical normal form 
is effectively a function of d integrals, see Equation (j2.83p . 

In the previous section we saw that in the monomials that form the symbol of a 
Hamilton operator in A^*'^ order quantum normal form the coordinate pairs {xk,Ck) 
(or equivalently {qk,Pk)), k = 1, . . . ,d, again have equal integer exponents. Hence, 
the symbol is effectively a function of / = piqi, Jk = ^{Pk + Q'k)^ k = 2, . . . ,d. We will 
now show that analogously to (12.830 the Hamilton operator in A^**^ order quantum 
normal form is a function of the d operators 

/ := Op[/] , Jk := Op[Jfc] , k = 2,...,d, (3.90) 

see Equations (j3.15p and (j3.16p . To this end recall that the Hamilton operator in 
quantum normal form is localised near the equilibrium point. We will say that two 
operators Op[^] and Op[i3] are equal near a point z = {q,p) in phase space if their 
symbols A and B are equal in a neighbourhood of z. To indicate this we write 

Op[^] Op[B] . (3.91) 



Theorem 3. Let Op[HqJp\ be a Hamilton operator in N order quantum normal 
form with respect to an equilibrium point of its principal symbol of saddle- centre- ■ ■ - 
centre type, and assume furthermore that the frequencies • • • ji^d associated with 
the d — 1 centres are linearly independent over 'Z. Then there exists a polynomial 
K^Qj^p -.W^ of order [N/2] such that 



Op[Hqnf] =0 K^^Npi^, J2,---,Jd) 



(3.92) 



In this theorem denotes the integer part of N/2. The proof of this Theorem 

is based on the fohowing 

Lemma 14. Let I = pq, J = ^{p^ + q^), and I = Op[/], J = Op[J], respectively, 
then there are integers r„ ^ such that for any n G M, 

fc=0 ^ ^ 

and 

//ere [n/2] denotes the integer part of n/2, and the coefficients T^^k o.re determined 
by the recursion relation 

Tn+i,k = r„,fc + n^r„_i,fc_i fork>l (3.95) 

and Tnfl = 1- 

Proof. We start by considering the case of I = pq. The strategy will be to use the 
Weyl calculus to determine the symbol of Op [I"] as a function of /, and then to 
invert this relation. The symbol of Op[/"'] is := I * I * ■ ■ ■ * I , the n-fold star 
product of /. Using I = pq and the definition of the star product in p. 171) we find 
the recursion relation ^ 

I*r = + Q j n^I"-! . (3.96) 

This can be rewritten as 

Op[/"+i] = / Op[/"] - 0) Op[/"-i] , (3.97) 

which can be used to determine the := Op[/"] recursively. If we insert the ansatz 
(I3.93P into the recursion relation ()3.97p we find the recursion for the coefficients 

In order to show the validity of Equation ()3.94p we apply the same strategy and 
find instead of (j3.97p 

Op[ J"+i] = J Op[ J"] + 0) \^ Op[ J"-i] , (3.98) 

and inserting now (|3.94p as an ansatz into this equation leads again to the relation 
(j3.95p for the coefficients. □ 

We note that the closed formulae for /" and given in [Cre90] are not correct. 
We now prove Theorem EJ 

Proof of Theorem El . It follows from our construction that the symbol of a Hamilton 
operator in quantum normal form is near {q,p) = (0,0) a polynomial in I, J^, 
k = 2, . . . ,d that can be written in the following form: 



^SF = E^'^''^"^^''^2""---^r> (3.99) 
1=1 



where 2j/+2 X]a:=i "k;i < N, or equivalently ji+^k=i "k;i < N/2, for all / = 1, . . . , L. 
For Op[ffq^p] we thus find 

L 

Op[^Sf] = E 0P[^"^1 Op[J2°^1 • • • Op[J°^^'] . (3.100) 
If we insert the expansions from Lemma [H] into (j3.100p we obtain 

L [ai;i/2] [a2;i/2] K,/2] 



i=l fci=0 fc2=0 A:d=0 



2(A,i+...+fc. 



X j jai;;-2A:i ja2;i-2fc2 _ _ _ jad;!" 

(3.101) 

Since + Ylt=i ak;i < N/2 for alU = 1, . . . , L it follows that the RHS of (ICTTT) 
is a polynomial of order in I, J2, . . . , Jd- This polynomial defines the function 

We note that for 7i — > the polynomial ^Q^p tends to the polynomial i^^NF 
defined in p. 830 that gives the N^^ order classical normal form as a function of the 
integrals / and Jf., k = 2, . . . ,d. Though this is obvious from the proof of Theorem [3] 
it is worth mentioning that in general the coefficients in the polynomial K^^Jp differ 

from the polynomial that is obtained from writing as a function of / and Jfc, 

k = 2, . . . , d. We will see this in the example presented in Sec. 13.51 

Theorem [3] is a crucial result. It tells us that the truncated quantum normal 
form simply is a polynomial in the operators / and J^, k = 2, . . . ,d, whose spectral 
properties are well known. As we will see in more detail in Sections [5] and [6] this 
will allow us to compute quantum reaction rates and quantum resonances with high 
efficiency. 



3.5 Quantum normal form for one-dimensional potential 
barriers 

In the following we present the explicit computation of the quantum normal form for 
Hamilton operators of one-dimensional systems of type "kinetic plus potential" where 
the potential has a maximum. It is important to point out that the applicability of 
the normal form algorithms - both classical and quantum - are not restricted to 
systems of the form "kinetic plus potential" (i.e., for example Coriolis terms in the 
Hamiltonian function or Hamilton operator due to a magnetic field or a rotating 
coordinate frame are allowed). Since even for this simple one-dimensional problem 
the expressions for the symbols and operators involved soon become very lengthy we 
will carry out the quantum normal form algorithm only to order 4. We note that 
we implemented the normalisation algorithm in the programming language C-|— |-. 
In our object oriented implementation the number of dimensions and the order of 
truncation of the normal form can be chosen arbitrarily. This C-|--|- program will 
be used to compute the high order quantum normal forms for the more complicated 
examples given in Section [71 

For now let us consider a Hamilton operator of the form 



where the potential V is assumed to have a (non-degenerate) maximum at q = qq. 
The Weyl symbol of H is given by 



H{h,q,p) = ^p^ + V{q), (3.103) 
Zni 

i.e., Op[i/] = H. Since the symbol H does not depend on h, the symbol agrees with 
the principal symbol. Hamilton's equations for the Hamiltonian function given by H 
then have an equilibrum point at {q,p) = {qo, 0) which is of saddle stability type, i.e., 
the matrix associated with the linearisation of the Hamiltonian vector field about 
the equilibrium point has a pair of real eigenvalues itA. Here A is given by 



X = J-lv"iqo). (3.104) 
V m 

The first two steps in the sequence of transformations ()3.65p serve to shift the 
equilibrium point of the (principal) symbol to the origin and to simplify the quadratic 
part of the symbol. As mentioned in Sec. 13.31 it follows from Lemma [6] (exact Egorov) 
that the transformations of the symbol H to achieve these goals agree with the 
corresponding classical transformations. 

Classically, we shift the equilibrium point to the origin of the coordinate system 
by transforming the coordinates according to 

{q,p) ^ {q - qo,p) ■ (3.105) 

For completeness, we note that this transformation can be obtained from the time 
one map of the flow generated by the first order polynomial 

Wi{q,p)=-qop, (3.106) 

i.e., ^lY^{q,p) = {q — qo,p)- The Weyl quantisation of Wi is given by 

Op[Wi] = qoih^. (3.107) 
aq 

It follows from Lemma that 
has the symbol 

{h, q,p)=Ho cD^i [h, q,p) = H{h, q + qo,p) = ^p" + V{q + go) • (3.109) 

We now want to find a unitary transformation such that the quadratic part of 
the symbol if*^^) of the transformed Hamilton operator assumes the form 

Hf\h,q,p) = \pq. (3.110) 
Classically, this is achieved by the transformation 

{q,p) ^ [y/m\q,—^p\ (3.111) 



'rnX 



followed by the 45° rotation 



{q,p)^ [l={p + q)^l^{p-q)\. (3.112) 



Both these transformations are symplectic. 

Again for completeness, we note that the transformation (j3.11ip can be obtained 
from the time one map of the flow generated by 

W2{q,p) = In {V^)pq, (3.113) 

i.e., 

'^w,{Q,p) = (V^q, -^p) ■ (3.114) 
\ vmA / 

The Weyl quantisation of W2 is given by 

Op[W2] = In {V^) J (^q^ + ^) , (3.115) 

see Equation (j3.16p . The transformation ()3.112p can be obtained from the time one 
map of the flow generated by 

W^(q,p) = j^{q^+p^), (3.116) 

which gives 

^Wiiq,P) =(^-^{p + q), ^{p - q)^ , (3.117) 

see the example after Lemma (exact Egorov), (j3.33p . The Weyl quantisation of 
is given by 

see Equation (j3.15p . 

Using Lemma [6] it follows that the symbol of 

Op[if(2)] = ei OPI^^I ei opt^^l Op[//(i)] e'i OpI^^I e'i OpI'^^] (3.119) 

is given by 

(2) (^n, q, p) = //(I) o o cD^l {h, q, p) 

00 k 00 

= Vo + Xqp+Y,Yl Vn^k^nP^'q'-" =■■ Yl ^kh^ I^P) > 

fc=3 n=0 fc=0 

where 

H^i\h,q,p)=Vo:=V{q^), H[^\h, q,p) = , H^^\h,q,p) = Xpq . (3.121) 
The coefficients of the monomials in (|3.12Up of cubic or higher degree are 

y . - (_i)n J \ d"^^'^(go) , + , > 3 (3 122) 

'^"'^ " ^ ^> n\j\ (2mA)("+^)/2 dq^+^ ' " + ^ ^ ^ . (3.122j 

So far, i.e., up to order 2, the transformations involved in the quantum normal 
form algorithm agree with their counterparts in the classical normal form algorithm. 
We now want to study the next steps in the sequence (j3.65p which give the quantum 
normal form of order three and four. To make these transformations well defined 
we from now on assume that we use the scheme outlined in Sec. l3.1.3l to localise the 
Hamilton operator H^'^^ and the operators which will generate the required unitary 



(3.120) 



transformations about the origin. The monomials in the third and fom'th order 
polynomials -^3^^ and H^'' have coefficients 

V,, = -y,3 = -lr,, = iv,, = -i^^y"'to), (3,123) 



respectively, where the primes denote derivatives. 

It follows from Equation (j3.70p that for W-^ G ^qm ioc symbol of the trans- 
formed operator 

Op[//(3)] = 0P[^3] Op[/f(2)]e-i Op[^3] (3.125) 

is given by 

^(3) ^ ^^3) ^ ^(3) ^ ^^3) ^ ^^3) ^ ^^3) ^ _ _ ^ ^3_^26) 

where following Lemma W2[ the terms Hj, and //^ agree for k < 2, and 

= + MadvFa = + {^^3, ^} , (3.127) 

//f ) = //f ) + Madiyg //f ^ + ^ [Madv^g]^ i^f ^ . (3.128) 

Equation (I3.127P is the homological equation. Introducing the operator 

V = {H^^\-} (3.129) 

the homological equation takes the form 

hP = hP-VW3, (3.130) 

which agrees with the form of the homological equation in Lemma [131 Following 
Sec. p.4p we need to solve the homological equation, i.e., choose W3, such that 
Since Wq^i = I™ ^|w3 ^ equivalently Ker Vl^^^ = {0}, we have to 

/ON 

choose W3 such that H^ ' = 0. Prom (fXT^fp we see that this is achieved by setting 
Ws{h,q,p) = -E A(2n-3) ^"-'-"^"^'""' ^^"^^^^ 

n=0 

= (^^' - - W + g') • (3.132) 

Inserting this W3 into ()3.128p gives 

y2 (3.133) 
(3/ + 12/g - 30/g2 + 12pg3 + 3g^ - Ah^) . 

A 

Note the occurrence of the term inolving h?. It is a consequence of the second term 
on the right hand side of (|3.128p which involves the Moyal bracket of two polynomials 
which are of degree higher than two for which the Moyal bracket no longer coincides 
with the Poisson bracket. 

Using Equation (|3.7Up again we see that for W4 S VV^j^. j^^,, the symbol of the 
transformed operator 

Op[i/(^)] = ei Op[l^4] Op[/f(3)]e-i Op[l4^4] (3.134) 



is given by 

= ^(4) + ^(4) + ^(4) ^ ^(4) + ^(4) ^ ^ _ _ ^ (3^^35) 

where it again follows from Lemma [T^ that H^"^ = H^^^ for A; < 3. For A; = 4 we 
obtain the homological equation 

//f ) = //f ) + Madvi/4 H^^^ = - ^^^4 ■ (3.136) 

We need to choose PV4 such that VH^^^ = 0. We therefore decompose H^^ according 
to 

= <L + <im> (3.137) 

/o\ 1 /ON I 

where -f^4.Ker ^ ^ w4 • -'-^ follows from Sec. 13.41 that 

(3) (3) 

-^4-Ker consists of all monomials of in which p and q have the same integer 

(3) (3) 

exponent and H^.^^^ consists of all monomials of in which p and q have different 
integer exponents. We thus have 

and 

^ifim = ^4;0 (p' - ^P^q - W + - ^ (3/ + 12^3^ + 12pg3 + 3/) . (3.139) 



To achieve Pi^f ^ = we choose 



VsjO /„ 4 , o/l^3„ o/l„„3 Q„4N ^4;0 / 4 n 3„ , q„„3 i\ 



Wi{h, q,P) = ^ + 24/<7 - 24pg^ - 3g4) - ^ (p^ _ Sp^g + 
We thus get 



(3.140) 



y2 

H^^\h, q,p) =Vo + Xpq + 6 V^^Qp^q^ + ^ (30/?^ + 4;,2) ^ ^ ^g^^^^^ 

A 

where the remainder O5 is defined according to Definition HI Neglecting O5 gives 
the symbol of the 4*'^ order quantum normal form. In order to get the corresponding 
operator we have to replace the factors I = pq hy the operator I = Op[/]. To this 
end we use the recurrence ()3.93p in Lemma [H] to get 

Op[/2] =P-^. (3.142) 

The 4**^ order quantum normal form of the operator H in ()3.102p is thus given by 

^qnf(^) = ^o + Ai+ I 30^+6y4;oJ Y f 7^ + 3y4;oJ (3.143) 

1 f 7 



(^(^'"te))^ + ^""(<Zo))/.^ (3.144) 



This gives the first correction term to the well known quadratic approximation which 
consists of approximating the potential barrier by an inverted parabola. The corre- 
sponding classical normal form is given by 



K^Un =Vo + XI + [^(^'"M + y"'\^o) ) I' . (3.145) 



We see that the polynomials i^q^p (in /) has two more terms than the polynomial 

Kq^p (in /) . These are the terms involving h, and their occurrence is due to the Moyal 
bracket (see the comment following Equation (j3.133p ) which enters the quantum 
normal form computation on the level of the symbols and the Weyl quantisation of 
powers of the classical integral I = pq (see ()3.142p ) which is required to obtain the 
Hamilton operator from its symbol. 



4 Classical Reaction Dynamics and Reaction 
Probabilities 

In this section we give an overview of the theory of reaction dynamics that is firmly 
rooted in the dynamical arena of phase space and has recently been developed in 
IWWJUOll IUJP+01[ IWBW04bl IWWOll IWBWOSal IWBWOSbl IWBWOScj . This sec- 



tion is organised as follows. In Sec. 14.11 we describe the geometric structures in 
phase space near an equilibrium point of saddle-centre-- • • -centre stability type (see 
Sec. 12. 3|) that control the classical dynamics of reactions. These phase space struc- 
tures are "realised" through the classical normal form, and details of this are given in 
Sec. 14.21 where we also provide a detailed discussion of how these phase space struc- 
tures constrain trajectories of Hamilton's equations. In Sec. 14.31 we describe how the 
integrability of the truncated normal form gives rise to the foliation of the phase 
space near the saddle by Lagrangian manifolds. This will be of central importance 
for the quantum mechanics of reactions as we will see in Sec. [5l In Sec. 14.41 we will 
show how the normal form can be used to compute the directional flux through the 
dividing surface. As we will see the normal form obtained from truncating the nor- 
mal form algorithm at a suitable order gives a very accurate description of the local 
dynamics. Means to verify the accuracy are discussed in Sec. 14.51 While the normal 
form technique is "locally applicable" in a neighbourhood of the reaction region, in 
Sec. 14.61 we discuss how the local structures mentioned above can be globalised in a 
way that their influence on reactions outside this "local" region can be determined. 
Finally, in Sec. 14.71 we comment on the flux-flux autocorrelation function formalism 
to compute classical reaction probabilities that is frequently utilised in the chemistry 
literature, its relation to our phase space theory, and the computational benefits of 
our approach over the flux-flux autocorrelation function formalism. 



4.1 Phase Space Structures that Control Classical Re- 
action Dynamics: An Overview of the Geometry 

Our starting point is an equilibrium point of Hamilton's equations of saddle-centre- 
• • • -centre stability type. Near (and we will discuss what we mean by "near" in 
Section 14. 5p such equilibrium points there exist lower dimensional manifolds that 
completely dictate the dynamics of the evolution of trajectories from reactants to 
products (or vice- versa). The normal form theory developed in Sec. [2] provides a 
transformation to a new set of coordinates, referred to as the normal form coordinates, 
in which these manifolds can be identified and explicitly computed, and then mapped 
back into the original, "physical" coordinates via the normal form transformation. 
In this section we give a brief description of these phase space structures, and in 
Sec. 14.21 we will describe how they constrain trajectories. 

We let Eq denote the energy of the saddle, and we consider a fixed energy E > Eq 
(and "sufficiently close" to Eq). We will also restrict our attention to a certain 



neighbourhood U, local to the equilibrium point. We will defer a discussion exactly 
how this region is chosen to Sec. 14.51 suffice it to say for now that the region is chosen 
so that an integrable nonlinear approximation to the dynamics yields structures to 
within a given desired accuracy. 

Near this equilibrium point the {2d — l)-dimensional energy surface in the 2(i- 
dimensional phase space K,^'^ has the structure of a "spherical cylinder" S'^'^^'^ x M, 
i.e., the Cartesian product of a {2d — 2) -dimensional sphere S'^'^~'^ and a line IR. 
The dividing surface that we construct locally separates the energy surface into two 
components; "reactants" and "products". This dividing surface which we denote by 
S'^~'^{E) has the structure of a {2d — 2)-dimensional sphere S"^^'"^. It can be shown 
to have the following properties: 

• The only way that trajectories can evolve from the reactants component to the 
products component (and vice- versa), without leaving the local region U, is by 
crossmg S]^~'^{E). We refer to this property of Sl'^~'^{E) as the "bottleneck 
property" 

• The dividing surface that we construct is free of local recrossings; any trajectory 
which crosses Slt~'^{E) must leave the neighbourhood U before it might possibly 
cross Sl'^~'^{E) again. 

• A consequence of the previous property of the dividing surface is that it mini- 
mizes the (directional) flux. It is thus the optimal dividing surface sought for 
in variational transition state theory [WW04j . 

The dividing surface S'^f~'^{E) itself is divided into two hemispheres: the forward 
reactive hemisphere B'^~^^{E), and the backward reactive hemisphere B'^^~^{E). The 
hemispheres B'^'^~f'^{E) and B'^~^{E) are topological (2(i — 2)-balls. These two hemi- 
spheres are separated by the equator of S'^~'^{E), which itself is a sphere of dimension 
{2d - 3). On Blf^^^{E) and bI^~^{E) the Hamiltonian vector field is transversal to 
each of these surfaces. This transversality is the mathematical manifestation of "no 
recrossing". Heuristically, "transversal" means that the Hamiltonian vector field 
"pierces" the surfaces, i.e., there is not point where it is tangential to the surface. 
Now the Hamiltonian vector field pierces the surfaces B'^~^^{E) and B'^~^{E) in 
opposite directions. Since the vector field varies smoothly from point to point, it 
must be tangential to the equator on which B'^~^^{E) and B^'^~^{E) are joined. 
More mathematically, the fact that the Hamiltonian vector field is tangential to the 
equator means that the equator is an invariant manifold. In fact, it is a so-called nor- 
mally hyperbolic invariant manifold (NHIM) |Wig94| , denoted by S^^il^{E), where 
normal hyperbolicity means that the expansion and contraction rates transverse to 
the manifold dominate those tangent to the manifold, and there are an equal number 
of independent expanding and contracting directions transverse to the manifold at 
each point on the manifold. This implies that it is "saddle like" in terms of stability 
(in our set-up there is one expanding direction and one contracting direction normal 
to the NHIM at each point on the NHIM). Heuristically, one can think of it as a "big 
saddle like surface". In fact, the {2d — 3)-dimensional NHIM is the energy surface of 
an invariant subsystem which has d—1 degrees of freedom, i.e., one degree of freedom 
less than the full system. In chemistry terminology this subsystem is the "activated 



^ Here we inserted the restriction 'without leaving the local region to exclude the case where the 
dividing surface does not divide the full (global) energy surface into two disjoint components. For example, 
two regions in an energy surface might be connected by channels associated with two different saddle-centre- 
• • • -centre equilibrium points. 



complex" , which may be thought of as representing an oscillating (unstable) "super- 
molecule" poised between reactants and products |Eyr35 , IPecSH IMil98aj . 



Normally hyperbolic invariant manifolds have stable and unstable manifolds, 
which themselves are invariant manifolds. In particular, the NHIM, 'S'^him(-^)' 
{2d — 2) -dimensional stable and unstable manifolds W'^{E) and W'^{E) which are 
isoenergetic, i.e., contained in the energy surface. These invariant manifolds have 
the topology of spherical cylinders 5^*^"^ x IR. Since they are of codimension one in 
the energy surface, i.e., they are of one dimension less than the energy surface, they 
act as impenetrable barriers. The importance of these particular geometrical struc- 
tures is that all reactive trajectories (both forward and backward) must lie inside 
regions of the energy surface that are enclosed by the NHIM's stable and unsta- 
ble manifolds. This can be described more precisely by first noting that W'^{E) and 
W^{E) each have two branches that "join' at the NHIM. We call these branches of the 
forward and backward branches of W'^{E) and H^"(£'), and denote them by Wj{E), 
W§{E), WJ{E) and W^{E), respectively. We call the union of the forward branches, 
Wf{E) := Wj{E)UWj{E), the forward reactive spherical cylinder. Trajectories with 
initial conditions enclosed by Wf{E) in the reactants component of the energy sur- 
face evolve towards the forward hemisphere of the dividing surface B'^~^^{E), cross 

-^ds7^(-^)' evolve into a region of the products component of the energy surface 
that is enclosed by Wf{E). Similarly, we call the union of the backward branches, 
Wi){E) := W^{E) U ly^(S), the backward reactive spherical cylinder. Trajectories 
with initial conditions enclosed by Wi,{E) in the products component of the energy 
surface evolve towards B'^'^~^{E), cross B'^~^{E), and evolve into a region of the 
reactants component of the energy surface that is enclosed by Wii{E). All forward 
reactive trajectories are enclosed by Wf{E) and all all backward reactive trajectories 
are enclosed by Wi,{E). As we will see in the next section these structures can be 
computed from the normal form developed in Section [2j 

4.2 The Normal Form Coordinates: Phase Space Struc- 
tures and Trajectories of Hamilton's Equations 

We now describe how the phase space structures mentioned in the previous section 
can be identified and computed from the normal form algorithm, and how they influ- 
ence trajectories of Hamilton's equations. From the discussion in Section 12.31 after 
steps of the normal form algorithm, we have constructed a coordinate transfor- 
mation from the original, "physical" coordinates to new, "normal form" coordinates 
(gj;^'' , . . . , q~^'^ , p^^^ , • • • , P^^'^ ) , and in these new coordinates the Hamiltonian trun- 
cated at order A'^ takes the form 

it(A^) _ j^W (riN) j{N) AN). 

-"CNF ~ -^CNFV-' ^■''2 ^---^-^d ) -j^-j 

= Eo + A/(^) + W2 Ja^-* + . . . + L^djji^^ + higher order terms , 

where 

/'"'=rpS"'. 4"'4((«r)'+(fn')- ^=2.....<i. (4.2) 

and the higher order terms are of order greater than 1 and less than or equal to [N/2] 
(in the integrals), see Equation (12.830 in Sec. 12.31 

The quantities (|4.2|) are integrals of the motion ("conserved quantities"), i.e, they 
are constant on trajectories of the Hamiltonian vector field given by the A^*'^ order 



classical normal form Hamiltonian Henceforth we will drop the superscripts (N) 
for the sake of a less cumbersome notation, but it should be understood that the 
normal form procedure is truncated at some fixed order N. 

In the normal form coordinates, using (|4.1|) and (|4.2|) . Hamilton's equations take 
the form 



Pl= -^^(/,J2,. 

dKcNF I T T„ 



, Jd) qi = 


A(/,J2,.. 


■ ■,Jd) qi 


,Jd)pi = 


-A(/,J2,.. 


■ , Jd) pi 


, Jd) Pk = 


^k{I, J2, ■ ■ 


■ , Jd)pk 


, Jd) Qk = 


-^kil, J2, ■ ■ 


■ ,Jd) qk 



(4.3) 



qk = 

Pk= -^^{I,J2,---,Jd)qk^ -Qk{I,J2,...,Jd)qk, k = 2,...,d. 



These equations appear "decoupled" . It is important to understand this statement in 
quotations since the equations are not "decoupled" in the usual fashion. Nevertheless, 
effectively, this is the case since the coefficient A and the nonlinear frequencies O^, 
k = 2, . . . ,d, are constant on a given trajectory. This follows from the fact that 
they are functions of the integrals / and Jk, k = 2,. . . ,d. Hence, once the initial 
conditions for a trajectory are chosen, then the coefficients of (|4.3p are constant (in 
time), and in this sense the equations are decoupled and can be easily integrated. 
The reason we have this property is a result of the d independent integrals given 
in (14. 2|) . However, A and the nonlinear frequencies Qk, k = 2,. . . ,d, will generally 
vary from trajectory to trajectory and the equations are hence not decoupled in the 
classical sense. We could view them as being "decoupled on trajectories" as a result 
of the d integrals being constant on trajectories. In mathematical terms this means 
that the equations of motion are integrable. The notion 'integrability' can be viewed 
as a generalisation of the notion 'separability'. The latter refers to the property of 
the equations of motion that allows to achieve a decoupling of the form (j4.3p from 
a transformation that involves the configuration space variables q only (which then 
entails a transformation of the momenta p to give a symplectic transformation of the 
full phase space coordinates). Historically, separability has played an important role 
in developing approximate transition state theory and analyzing tunnelling effects, 
see, e.g., |JR6H IEM741 IMil76| . Indeed, if the full dynamics is separable near the 
saddle point (in phase space) then the construction of a dividing surface with no 
recrossing is trivial and the choice of reaction coordinate is "obvious". However, it 
is important to point in the neighbourhood of a saddle-centre-- • • -centre equilibrium 
point the equations of motions are in general not separable but the normal form 
transformation leading to the decoupling in ()4.3p in general involves a symplectic 
transformation which mixes configuration and momentum variables. 
The general solution of (14. 3p is given by 



gi(t) = Aiexp( A(I, J2,..., Jd)t) , 

pi{t) = Bi exp ( - A(I, J2, . . . , Jd) t) , 

Qkit) = Ak sin {Qkil, J2, ■ ■ ■ ,Jd)t + (pk) , 

Pk{t) = Ak cos {Qk{I, J2,---,Jd)t + (pk), /c = 2, . . . , d , 

where the Ai, . . . , Ad, ip2, ■ ■ ■ ,^Pd and Bi are 2d constants determined by the initial 
conditions (gi(0), . . . , grf(0),pi(0), . . . ,pd(0)). The constants in (j4.4p determine the 



^Tlie fact that there are d constants of motion is a consequence of the non-resonance assumption on 
the linear frequencies Uk, k = 2, . . . , d. If there are resonances amongst the cok, then there will be fewer 
integrals. 




Figure 2: Projection of energy surfaces (turquoise regions) to the planes of the normal form coordinates. 
The energy surface in the top panel has E < Eq; the energy surface in the bottom panel has E > Eq- 



integrals according to 

I = AiBu Jk = \Al, k = 2,...,d. (4.5) 

From the general solution ()4.4p we see that the motion is generally hyperbolic (i.e., 
"saddle like") in the plane of the coordinates {qi,pi) associated with the saddle and 
rotational in the planes of the coordinate pairs {qk,Pk), k = 2, . . . , d, associated with 
the centre directions. 

In the following, we show how the normal form, which is valid in the neighbour- 
hood of the saddle-centre- • • -centre equilibrium point, gives explicit formulae for the 
various manifolds described in Sec. 14.11 At the same time, we show how trajectories 
of Hamilton's equations expressed in the normal form coordinates, are constrained 
by these manifolds. Many more details can be found in [U.TP+nH IWBWn4bj . The 
geometrical illustrations that we give are for three degrees of freedom. In fact, con- 
ceptually, the step from two to three degrees of freedom is the big step; once the 
case of three degrees of freedom is well understood, it is not difficult to incorporate 
more degrees of freedom. We begin by describing the local structure of the energy 
surfaces. 

The structure of an energy surface near a saddle point: For E < Eq, 

the energy surface consists of two disjoint components. The two components cor- 
respond to "reactants" and "products." The top panel of Fig. [2] shows how the 
two components project to the various planes of the normal form coordinates. The 
projection to the plane of the saddle coordinates {qi,pi) is bounded away from the 
origin by the two branches of the hyperbola, qipi = I < 0, where / is given implicitly 
by the energy equation with the centre actions Jk, k = 2, . . . ,d, set equal to zero: 
-f^CNF(-^, 0, . . . ,0) = E < Eq. The projections to the planes of the centre coordinates, 
{qk,Pk), k = 2, . . . ,d, are unbounded. 

At E = Eq, the formerly disconnected components merge (the energy surface bi- 
furcates), and for E > Eq the energy surface has locally the structure of a spherical 
cylinder, S'^'^^^ x El. Its projection to the plane of the saddle coordinates now in- 
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Figure 3: Projection of the dividing surface and reacting and nonreacting trajectories to the planes 
of the normal form coordinates. In the plane of the saddle coordinates, the projection of the dividing 
surface is the dark red diagonal line segment, which has q\ = pi- In the planes of the centre coordinates, 
the projections of the dividing surface are the dark red discs. Forward and backward reactive trajectories 
(yellow and blue) project to the first and third quadrant in the plane of the saddle coordinates, respectively, 
and pass through the dividing surface. The red and green curves mark nonreactive trajectories on the 
reactant side {pi — qi > 0), and on the product side (pi — gi < 0), of the dividing surface, respectively. 
The turquoise regions indicate the projections of the energy surface. 



eludes the origin. In the first and third quadrants it is bounded by the two branches 
of the hyperbola, qipi = / > 0, where / is again given implicitly by the energy equa- 
tion with all centre actions equal to zero, but now with an energy greater than Eq: 
KcNpil, 0, . . . ,0) = E > Eq. The projections to the planes of the centre coordinates 
are again unbounded. This is illustrated in the bottom panel of Fig. [2l 

The dividing surface and reacting and nonreacting trajectories: On an 

energy surface with E > Eq, we define the dividing surface by qi = pi- This gives a 
(2(i— 2)-sphere which we denote by S'^^~'^{E). Its projection to the saddle coordinates 
simply gives a line segment through the origin which joins the boundaries of the 
projection of the energy surface, as shown in Fig. [3l The projections of the dividing 
surface to the planes of the centre coordinates are bounded by circles (p|-|-g^)/2 = J^, 
k = 2, . . . ,d, where Jk is determined by the energy equation with the other centre 
actions, J/, I ^ k, and the saddle integral, /, set equal to zero. The dividing surface 
divides the energy surface into two halves, Pi — qi > and Pi—qi < 0, corresponding 
to reactants and products. 

As mentioned above, trajectories project to hyperbolae in the plane of the saddle 
coordinates, and to circles in the planes of the centre coordinates. The sign of / 
determines whether a trajectory is nonreacting or reacting, see Fig. [3l Trajectories 
which have / < are nonreactive and for one branch of the hyperbola qipi = I they 
stay on the reactants side and for the other branch they stay on the products side; 
trajectories with / > are reactive, and for one branch of the hyperbola qipi = I 
they react in the forward direction, i.e., from reactants to products, and for the other 
branch they react in the backward direction, i.e., from products to reactants. The 
projections of reactive trajectories to the planes of the centre coordinates are always 
contained in the projections of the dividing surface. In this, and other ways, the 
geometry of the reaction is highly constrained. There is no analogous restriction on 
the projections of nonreactive trajectories to the centre coordinates. 

The normally hyperbolic invariant manifold (NHIM) and its relation to 
the 'activated complex': On an energy surface with E > Eq, the NHIM is given 
by qi = pi = 0. The NHIM has the structure of a (2d — 3)-sphere, which we denote 
by 5'NfflM(^)- The NHIM is the equator of the dividing surface; it divides it into 



two "hemispheres": the forward dividing surface, which has qi = pi > 0, and the 
backward dividing surface, which has qi = pi < 0. The forward and backward 
dividing surfaces have the structure of {2d — 2)-dimensional balls, which we denote 
by B'^~^^{E) and B'^~^{E), respectively. All forward reactive trajectories cross 

B'^~^^{E); all backward reactive trajectories cross B'^~^{E). Since gi = pi = in 
the equations of motion (j4.3p implies that = pi = 0, the NHIM is an invariant 
manifold, i.e., trajectories started in the NHIM stay in the NHIM for all time. The 
system resulting from gi = pi = is an invariant subsystem with one degree of 
freedom less than the full system. In fact, gi = pi = defines the centre manifold 
associated with the saddle-centre-- • • -centre equilibrium point, and the NHIM at an 
energy E greater than the energy of the quilibrium point is given by the intersection 
of the centre manifold with the energy surface of this energy E lu.TP+niilwWni] . 



This subsystem is the "activated complex" (in phase space), located between 
reactants and products (see Sec. 14. ip . The NHIM can be considered to be the energy 
surface of the activated complex. In particular, all trajectories in the NHIM have 
/ = 0. 

The equations of motion (14. 3|) also show that pi — qi < on the forward divid- 
ing surface B'^~^^{E), and pi — gi > on the backward dividing surface B'^'^~^{E). 
Hence, except for the NHIM, which is is an invariant manifold, the dividing surface is 
everywhere transverse to the Hamiltonian flow. This means that a trajectory, after 
having crossed the forward or backward dividing surface, B'^~^^{E) or B'^~^{E), 
respectively, must leave the neighbourhood of the dividing surface before it can pos- 
sibly cross it again. Indeed, such a trajectory must leave the local region in which 
the normal form is valid before it can possibly cross the dividing surface again. 

The NHIM has a special structure: due to the conservation of the centre actions, 
it is filled, or foliated, by invariant {d — l)-dimensional tori, T'^~^. More precisely, 
for d = 3 degrees of freedom, each value of J2 implicitly defines a value of J3 by the 
energy equation i^cNF(0, J2, ^3) = E. For three degrees of freedom, the NHIM is 
thus foliated by a one-parameter family of invariant 2-tori. The end points of the 
parameterization interval correspond to J2 = (implying q2 = P2 = 0) and J3 = 
(implying qs = Ps = 0), respectively. At the end points, the 2-tori thus degenerate 
to periodic orbits, the so-called Lyapunov periodic orbits. As we will discuss in more 
detail in Sections [5] and [U the fact that the NHIM is foliated by invariant tori has 
important consequences for the corresponding quantum system. 

The stable and unstable manifolds of the NHIM forming the phase space 
conduits for reactions: Since the NHIM is of saddle stability type, it has stable 
and unstable manifolds, W^{E) and W^{E). The stable and unstable manifolds have 
the structure of spherical cylinders, 5'^°'"^ x R. Each of them consists of two branches: 
the "forward branches", which we denote by W^{E) and W^{E), and the "backward 
branches", which we denote by W^{E) and W^{E). In terms of the normal form 
coordinates, W^{E) is given by gi = with pi > 0, Wj{E) is given by pi = with 
qi > 0, W^{E) is given by gi = with pi < 0, and W^{E) is given by = with 
qi < 0, see Fig. [H Trajectories on these manifolds have / = 0. 

Since the stable and unstable manifolds of the NHIM are of one less dimension 
than the energy surface, they enclose volumes of the energy surface. We call the union 
of the forward branches, Wj{E) and Wj{E), the forward reactive spherical cylinder 
and denote it by Wf{E). Similarly, we define the backward reactive spherical cylinder, 
Wb{E), as the union of the backward branches, W^{E) and W^{E). 

The reactive volumes enclosed by Wf{E) and Wh{E) are shown in Fig. as 
their projections to the normal form coordinate planes. In the plane of the saddle 




Figure 4: The projection of the NHIM and the local parts of its stable and unstable manifolds, W^{E) 
and W^{E), to the planes of the normal form coordinates. In the plane of the saddle coordinates, the 
projection of the NHIM is the origin marked by the blue bold point, and the projection of W''{E) and 
W"{E) are the pi-axis and gi-axis, respectively. W{E) consists of the forward and backward branches 
Wf{E) and wf{E), which have pi > and pi < 0, respectively; W^E) consists of Wf{E) and W^{E), 
which have gi > and qi < 0, respectively. In the plane of the centre coordinates, the projections of the 
NHIM, W''{E), and W^{E) (the blue circular discs) coincide with the projection of the dividing surface 
in Fig. [31 The turquoise regions mark the projections of the energy surface. 




Figure 5: Projections of the reactive volumes enclosed by the forward and backward reactive spherical 
cylinders, Wf{E) and Wii{E), and the forward and backward reactions paths, to the planes of the normal 
form coordinates. The volumes enclosed by Wf{E) and Wb{E) project to the dark pink and green regions 
in the first and third quadrant in the plane of the saddle coordinates, respectively. These volumes project 
to the dark green/dark pink brindled disks in the planes of the centre coordinates, where their projections 
coincide with the projection of the NHIM and the dividing surface in Figs. [3] and [H The forward and 
backward reaction paths project to the two branches of a hyperbola marked blue in the first and third 
quadrant in the plane of the saddle coordinates, respectively, and to the origins (bold blue points) in the 
planes of the centre coordinates. The turquoise regions mark the projections of the energy surface. 



coordinates, the reactive volume enclosed by Wf{E) projects to the first quadrant. 
This projection is bounded by the corresponding hyperbola qipi = /, with / obtained 
from Kcnf{I,0, . ■ ■ ,0) = E. Likewise, Wb{E) projects to the third quadrant in 
the (gi,pi)-plane. Wf{E) encloses a/Horward reactive trajectories; Wh{E) encloses 
all backward reactive trajectories. All nonreactive trajectories are contained in the 
complement. 

Forvi^ard and backward reaction paths: The local geometry of Wf{E) and Wh{E) 
suggests a natural definition of dynamical forward and backward reaction paths as 
the unique paths in phase space obtained by putting all of the energy of a reacting 
trajectory into the reacting mode, i.e., setting q2 = • • • = Qd = P2 = • • • = Pd = 
0. This gives the two branches of the hyperbola qipi = I, with / obtained from 
-f^CNF(-^, 0, . . . ,0) = E, which in phase space are contained in the plane of the saddle 
coordinates, see Fig. [5l This way, the forward (respectively, backward) reaction path 
can be thought of as the "centre curve" of the relevant volume enclosed by the forward 
(resp., backward) reactive spherical cylinder Wf{E) (resp., Wb{E)). These reaction 
paths are the special reactive trajectories which intersect the dividing surface at the 
"poles" (in the sense of North and South poles, where qi = pi assumes its maximum 
and minimum value on the dividing surface). 

The Transmission Time through the Transition State Region: The normal 
form coordinates provide a way of computing the time for all trajectories to cross the 
transition region. We illustrate this with a forward reacting trajectory (a similar ar- 
gument and calculation can be applied to backward reacting trajectories). We choose 
the boundary for the entrance to the reaction region to he pi — qi = c for some con- 
stant c > 0, i.e., initial conditions which lie on the reactant side of the transition state, 
and the boundary for exiting the reaction region to be Pi — qi = —c on the product 
side. We now compute the time of flight for a forward reacting trajectory with initial 
condition on pi — qi = c to reach Pi — qi = — c on the product side. The solutions 
are qi{t) = qi{0) exp(A(/, J2, . . . , Jd)t) and pi{t) = pi{0) exp(-A(/, J2, . . . , Jd)t) (see 
(|4.4p ). where A(/, J2, . . . , Jd) is determined by the initial conditions. This gives the 
time of flight as 



The time diverges logarithmically as qi{0) — > 0, i.e., the closer the trajectory starts 
to the boundary Wf{E). It is not difficult to see that the time of flight is shortest for 
the centre curve of the volume enclosed by Wf{E), i.e., the trajectory which traverses 
the transition state region fastest is precisely our forward reaction path. A similar 
construction applies to backward reactive trajectories. 

In fact, the normal form can be used to map trajectories through the transition 
state region, i.e. the phase space point at which a trajectory enters the transition state 
region can be mapped analytically to the phase space point at which the trajectories 
exits the transition state region. 

4.3 The Normal Form Coordinates: The FoUation of the 
Reaction Region by Lagrangian Submanifolds 

In Section 14.21 we have indicated that the different types of possible motion near a 
saddle-centre- • • -centre equilibrium point can be described in terms of the integrals. 
In fact, the existence of the d integrals ()4.2p leads to even further constraints on the 




(4.6) 



classical motions and hence to even more detailed structuring of the phase space near 
a saddle-centre-- • • -centre equilibrium point than we already described in Sec. 14. 2[ As 
we will see in Sect. [5] this structure will have important consequences for the quantum 
mechanics of reactions. In order to describe this structure we introduce the so called 
momentum map M. |Gui94t IMR99j which maps a point (gi, . . . , . . . ^pd) in the 

phase space K,^*^ to the integrals evaluated at this point: 

M{qi,. . . ,qd,Pi,. . . ,Pd) ^ {I,J2,---,Jd)- (4.7) 

The preimage of a value for the constants of motion (I, J2, . . . , Jd) under Ai is called 
a fibre. A fibre thus corresponds to the common level set of the integrals in phase 
space. 

A point {qi, . . . ,qd,Pi, ■ ■ ■ ,Pd) is called a regular point of the momentum map 
if the linearisation of the momentum map, DM, has rank d at this point, i.e., if 
the gradients of the d integrals I, Jk, k = 2, . . . ,d, with respect to the phase space 
coordinates {q,p) are linearly independent at this point. If the rank of DM is less 
than d then the point is called an irregular point. A regular fibre is a fibre which 
consists of regular points only. The regular fibres of the momentum map in ()4.7p are 
d-dimensional manifolds given by the Cartesian product of an hyperbola qipi = I in 
the saddle plane {qi,pi) and d—1 circles in the centre planes {qkiVk)-, k = 2, . . . , d. 
Since hyperbola qipi = I consists of two branches each of which have the topology of 
a line IR, the regular fibres consist of two disjoint toroidal cylinders, T*^"^ x IR, which 
are the Cartesian products of a (d— l)-dimensional torus and a line. We denote these 
toroidal cylinders by 

^?;j2,...,J, ={('?'^') elt''^ : Piqi=I, l{pl + ql) = J2,... ,\{pl + qd) =Jd,qi>0} 

(4.8) 

and 

^7,J2,...,J, = {('?'^') el^''^ ■ Pl1l=I^ l{pl + ql) =J2,... ,^{pl + qj) = Jd,qi <0}. 

(4.9) 

A^j^ Jd J2 Jd Lagrangian manifolds [Arn78j . Prominent examples 

of Lagrangian manifolds are tori which foliate the neighbourhood of a centre-. . .- 
centre equilibrium point and whose semiclassical quantisation often lead to a very 
good approximation of part of the energy spectra of the corresponding bounded 
system [OdABS]. In our case the Lagrangian manifolds are unbounded. They are the 
products of {d— 1) -dimensional tori T"^"^ and unbounded lines R. The toroidal base 
of these cylinders will again lead to semiclassical quantisation conditions and as we 
will see in Sections and E] this will have important consequences for the computation 
of quantum reaction rates and resonances. 

If the fibre contains an irregular point then the fibre is called singular. The image 
of the singular fibres under the momentum map is called the bifurcation diagram. It 
is easy to see that the bifurcation diagram consists of the set of (/, J2, . . . , Jd) where 
one or more of the integrals vanish. In Fig. [6]we show the image of the energy surface 
with energy E > Eq under the momentum map M in the space of the integrals for 
d = 3 degrees of freedom. The bifurcation diagram (of the energy surface) consists 
of the intersections of the image of the energy surface (the turquoise and green/dark 
ping brindled surface in Fig. E]) with one of the planes / = 0, J2 = or J3 = 0. 
Upon approaching one of the edges that have J2 = or J3 = the circle in the plane 
{Q2,P2) or ((73,^3), respectively, shrinks to a point, and accordingly the regular fibres 

X ]R reduce to cylinders or 'tubes' x IR,. At the top corner in Fig. [6] both J2 and 




Figure 6: Sketch of the image of the energy surface of energy E > Eq under the momentum map M. in 
Equation (14. 7p in the space of the integrals / and Jfc, = 2, . . . , d, for the case of d = 3 degrees of freedom. 
The green/dark pink brindled piece of the image of the energy surface has / > 0; the turquoise piece has 
/ < 0. The intersections with the planes / = 0, J2 = and J3 = (pieces of which are visualised by 
semitransparent planes for clarity) form the bifurcation diagram of the energy surface. The image of the 
energy surface is not bounded in the direction of negative I as indicated by the dashed line at the bottom. 
The topology of the fibres A^~^(/, J2, J3) is indicated for the various points {I,J2,J3) marked by dots. 
The fibre of a point (/, J2, J3) with / ^ consists of two disconnected manifolds as indicated by the factor 
of 2. The fibre of a point (/, J2, J3) with 1 = consists of a single connected manifold. 



J3 are zero. Here both circles in the centre planes ((72,^2) and (93,^3) have shrinked 
to points. The corresponding singular fibre consists of two lines, IR, which are the 
forward and backward reaction paths, respectively (see also Fig. [5]). 

The fibres mentioned so far all have / 7^ and each consist of a pair of two 
disconnected components. For / < 0, one member of each pair is located on the 
reactants side and the other on the products side of the dividing surface. For / > 0, 
one member of each pair consists of trajectories evolving from reactants to products 
and the other member consists of trajectories that evolve from products to reactants. 
In fact the two members of a fibre which has / > are contained in the energy surface 
volume enclosed by the forward and backward reactive spherical cylinders Wf{E) 
and Wb{E), see Fig. [5l For this reason we marked the piece of the image of the 
energy surface under the momentum map which has / > by the same green/dark 
pink colour in Fig. [6] that we used Fig. \5\ Green corresponds to forward reactive 
trajectories and dark pink corresponds to backward reactive trajectories. Under the 
momentum map these trajectories have the same image. 

The light blue line in Fig. [6] which has / = is the image of the NHIM under 
the momentum map. For three degrees of freedom the NHIM is a 3-dimensional 
sphere, and as mentioned in Sec. 14.21 and indicated in Fig. [6] it is foliated by a one- 
parameter family of invariant 2-tori which shrink to periodic orbits, i.e. circles S^, at 
the end points of the parameterisation interval. As we already indicated in Sec. 14.21 
this foliation of the NHIM has important consequences for the quantum mechanics of 
reactions which we will dicuss in Sections [S] and El Moreover, we will see in Sec. 14.41 
that, for d = 3 degrees of freedom, the area enclosed by the image of the NHIM in 
the plane {J2, J3) gives, up to a prefactor, the directional flux through the dividing 
surface. 




Figure 7: Contour i^cNF(0, J2, • ■ • , Jd) = E (blue line) in the space of the centre integrals (J2, . . . , Jd) 
for d = 3 degrees of freedom. Up to the prefactor {2tt)'^^^, the area V{E) of the enclosed region (marked 
green) gives the directional flux through the dividing surface, see Equation ()4.10p . The green region agrees 
with the projection of the piece of the image of the energy surface under the momentum map which has 
/ > in Fig. [6]to the (J2, J3)-plane. 



4.4 The Directional Flux Through the Dividing Surface 

A key ingredient of transition state theory and the classical reaction rate is the direc- 
tional flux through the dividing surface defined in Sec. 14.11 Given the Hamiltonian 
function in normal form expressed as a function of the integrals (|4.2|) . and a fixed 
energy E above the energy of the saddle-centre-- • • -centre, Eq, it is shown in |WW04| 
that the directional flux through the dividing surface is given by 

f{E) = {27ry-'V{E), (4.10) 

where V{E) is the volume in the space of the actions (J2, . . . , Jd) enclosed by the 
contour Hcnf{0, J2, ■ ■ ■ , Jd) = E. For E < Eq, the directional flux is zero. For the 
case of a system with three degrees of freedom for which we sketched the image of the 
energy surface in the space of the integrals in Fig. [6l the volume V{E) is given by the 
area in the (J2, J3) plane enclosed by the light blue line corresponding to the NHIM in 
Fig. m For clarity we illustrate this area again in Fig. [71 As we mentioned in Sec. 14.21 
the NHIM can be considered as the energy surface of an invariant subsystem with 
one degree of freedom less than the full system which is referred to as the activated 
complex in the chemistry literature. Therefore the flux can be interpreted as the 
volume enclosed by the energy surface (given by the NHIM) in the phase space of 
this invariant subsystem. This gives a direct connection between the directional flux 
through the dividing surface and the activated complex. In fact, the dimensionless 
quantity 

^Weyl(i^) = ^^{^, (4.11) 

where 2TTh is Planck's constant, is Weyl's approximation of the integrated density of 
states, or equivalently the mean number of quantum states of the activated complex 
with energies less than or equal to E (see, e.g. |Gut90j ). As we will see in Sec. [5] 
N'\j^Qyl{E) can be interpreted as the mean number of open quantum "transition chan- 
nels" at energy E. 

In the case where we only take into account the quadratic part of the normal form, 
or equivalently, if we linearise Hamilton's equations, we have HcNFil, J2, ■ ■ ■ , Jd) = 
XI + Ylk=2 ^kJk and the energy surface -ffcNF(0, J2, . . . , Jd) = encloses a simplex 



in (J2, . . . , Jd) whose volume leads to the well-known result [Mac90] 



This shows, e.g, that the flux scales with E'^~^ for energies close to the saddle energy. 
The key advantage of the normal form algorithm that we presented in Sec. [2] is that 
it allows one to include the non-linear corrections to (|4.12p to any desired order. 

Here we give a brief outline of the essential elements of the derivation of the 
expression for the flux in ()4.10p following the discussion |WW04j . It is important to 
note that our work is flrmly rooted in phase space. In particular, we are considering 
the (directional) flux of a vector fleld on phase space (Hamilton's equations) through 
a dividing surface in phase space (which has been proven to have the "no recrossing" 
property as discussed earlier). For this reason the modern notation of differential 
forms, especially in light of its importance in the modern formulation of Hamiltonian 
mechanics, proves to be most convenient and notationally economical. 

Therefore we begin by considering the phase space volume form Q = dpi A dqi A 
• • • A dpd A dqd, which in terms of the symplectic 2-from lu = Ylk=i '^Pk ^ '^Qk can be 
written as 17 = uj'^/dl. Note that in our case the phase space coordinates {q,p) used 
here will be A^**^ order normal form coordinates and we do not use superscripts {N) 
to indicate this. However the quantities introduced in the following do not depend 
on the chosen coordinate system. They are invariant under symplectic coordinate 
transformations. Let r] be an energy surface volume form defined via the property 
dH A = O. Then the flux through a codimension one submanifold of the {2d — 1)- 
dimensional energy surface H = E \s obtained from integrating over it the "flux" 
form Vt' given by the interior product of the Hamiltonian vector field Xh with ry 
|Mac90j . i.e. 

^' = ^^-'^=(7^1)1'^'"'' (4.13) 

where ixgiliCi, ■ ■ ■ ,6^- 2) = r]{^i, . . . ,^2d-2, Xh) for any 2d - 2 vectors Cfc- The 
second equality in (j4.13p is easily established on a non-critical energy surface, i.e. on 
an energy surface which contains no equilibria. The flux form Q' is exact. In fact the 
generalised "action" form 



d 

d-2 

(d-1)! 



2_^Pkdqk A , _ - 



k=l 

has the property d</> = 0' and facilitates the use of Stokes' theorem to compute the 
flux. In the case of two degrees of freedom we simply have Q' = uj = dpi A dqi + dp2 A 
dq2 and (p becomes the usual action form (p = pidqi+p2dq2 ■ Since the dividing surface 
S^^~'^{E) is a sphere, that is, a manifold without boundary, it follows from Stokes' 
theorem that the integral of O' over S'^f^^{E) is zero. In order to compute reaction 
rates one has to distinguish between the directions in which the Hamiltonian flow 
crosses the dividing surface (i.e., distinguish between forward and backward reactive 
trajectories). Given a normal bundle lj over Sl^-^iE) the direct ion can be specifled 
by the sign of the scalar product between the normal vectors and the Hamiltonian 
vector field. This scalar product is strictly positive on one hemispheres of S'^f^'^{E), 
strictly negative on the other hemisphere and zero only at the equator of S'^'^~'^{E), 



""^^Roughly speaking, at each point of the dividing surface we consider the normal vector in the energy 
surface. The normal bundle is the union of all vectors taken over all points on the dividing surface. 



i.e. at the normally hyperbolic invariant manifold S^^j'^{E), where the Hamiltonian 
vector field is tangent to S'^'^~'^{E). Likewise, the flux form 17' on S'^^~'^{E) vanishes 
nowhere on B'^~^^{E) and B'^~^{E) and is identically zero on S'^-^^{E). It is natural 
to take as the orientation of B'^~^{E) and B'^~^{E) the orientation they inherit 
from the dividing surface. Without restriction we may assume that the orientation 
of S'^~'^{E) is such that 0' is positive on the forward hemisphere B'^~^^{E) and 
negative on the backward hemisphere B'^'^~^{E), i.e. and —il.' can be considered 
as volume forms on B^~^{E) and B'^~^{E)^ respectively. It follows from Stokes' 
theorem that the flux through the forward and backward hemispheres, /j 



and 
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il.' , have the same magnitude but opposite sign and can be computed 



from integrating the action form (p over the NHIM: 
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(4.14) 



We call the positive quantity L 



'(E) 



ri' the forward flux and the negative quantity 



J„2d-2(^.Cl' the backward flux through S"^^ "^{E) ■ 

Writing the flux form 17' in terms of "angle-action variables" ((^i, . . . , ^pd^ I,J2,---, Jd) 
(these were derived in terms of the integrals of the normal form in Section 12.3. ip we 
obtain the result that the forward flux through the dividing surface is given by the 
expression in Equation (j4.10p . 



4.5 The Normal Form Coordinates: Issues Associated 
with Truncation 

The final question to address concerns the "validity" of the normal form transforma- 
tion. More precisely, this means how large can the neighbourhood (in phase space) 
U of the saddle-centre-- • • -centre equilibrium point be taken so that the geometric 
structures given by the normal form are accurate for the "full equations" . Actually, 
there are a number of questions to be answered related to "validity" . 

• In truncating the Taylor expansion of the Hamiltonian at degree A^, how do 
you determine A^? 

• What is the region of validity of the normal form transformation for the Taylor 
expanded Hamiltonian truncated at degree A^? 

• How "accurate" are the phase space structures (e.g., the dividing surface, the 
NHIM) for the normal form of the Hamiltonian truncated at degree A^? 

• How accurate are trajectories of the normal form of the Hamiltonian at order 

A^? 

First, general theory ensures us that the phase space structures exist, and have 
the properties described above (e.g., normal hyperbolicity, the bottleneck property, 
etc.), for energies sufficiently close to that of the saddle- centre- ■ ■ -centre equilibrium 
point |Wig94 . The normal form computation is merely an approach for realizing the 



geometrical structures that the theory tells us must exist. 

In practice, one Taylor expands the Hamiltonian and then truncates it at a degree 
that one thinks will provide sufficient accuracy for the range of energies of interest. 
Experience will generally provide some good "rules of thumb", e.g. for the HCN 
isomerization work described in |WBW04b] . an expansion up to degree 10 was found 



to provide sufficient accuracy in the range of energies studied (up to 0.2 eV above 
the saddle-centre-- • • -centre equiUbrium point). 

There is stih the question of accuracy. Once the normal form is computed to the 
desired degree (and, most importantly, the transformation and its inverse between 
the original coordinates and the "normal form coordinates"), and the energy is fixed, 
we have explicit formulae for the dividing surface, the NHIM (the "equator" of the 
dividing surface), and the (local) stable and unstable manifolds of the NHIlvf^. We 
next need to check their "accuracy". There are several tests that we employ, and 
these tests are carried out at fixed energy. 

• Numerically verify that the dividing surface satisfies the "bottleneck property" , 
i.e. it (locally) separates the energy surface into two components, and the only 
way a trajectory can pass between components (while remaining in this region) 
is by passing through the dividing surface. 

• Using the inverse of the normal form transformation map the NHIM and its 
(local) stable and unstable manifolds back into the original coordinates and 
check that the full (i.e., not a truncated Taylor expansion) Hamiltonian vector 
field is tangent to these surfaces. This is a requirement for these surface to be 
"invariant manifolds" . The tests are carried out pointwise on a grid of points 
covering the surfaces. 

• The integrals (14. 2p are constant in time on trajectories of the normal form of 
the truncated Taylor expansion. We check how they vary in time on trajectories 
of the full Hamiltonian. 

If the desired accuracy is obtained for this energy, then the energy may be in- 
creased and the accuracy tests are repeated at the higher energy. If accuracy is 
inadequate, then a higher degree Taylor expansion can be computed. As energy is 
increased, ultimately two factors may lead to break down of this approach for real- 
izing these phase space structures. One is that the energy surface may deform in 
such a way that the bottleneck property does not hold. Another is that the approach 
will require such a high degree Taylor expansion that it becomes computationally 
intractable. 

4.6 The Global Dynamics Associated with the Mani- 
folds Constructed in the Reaction Region 

As we have shown, the normal form transformation to normal form coordinates pro- 
vides a method for providing a complete understanding of the geometry of reaction 
dynamics in a neighbourhood U (in phase space) of the saddle-centre-- • • -centre equi- 
librium point of Hamilton's equations. By this, we mean that in the normal form 
coordinates we can give an explicit equation for the surfaces and, as a result of 
the "simple" structure of Hamilton's equations in the normal form coordinates, we 
can describe precisely the influence of these geometrical structures on trajectories of 
Hamilton's equations. In Tab. [I] we summarize the results obtained this far by pro- 
viding a list of the different surfaces that control the evolution of trajectories from 
reactants to products in the neighbourhood U in Fig. HI 

^^Here "local" means that we only have realizations of the stable and unstable manifolds in a neigh- 
bourhood of the saddle-centre-- • - -centre equilibrium point where the normal form transformation has the 
desired accuracy 



Geometrical Structure 


Equation in Normal Form Coordinates 


dividing surface, Sl^~'^{E) 


Qi = Pi 


forward reactive hemisphere, B'^~^^[E) 


gi = Pi > 


backward reactive hemisphere, B^1^{E) 


qi=pi <0 


NHIM, Sf^^{E) 


qi = Pi = 


stable manifold of the NHIM, W'{E) 


gi = 0, ^ 


unstable manifold of the NHIM, ^{E) 


= 0, gi ^ 


forward branch of ^{E), W^f{E) 


gi = 0, pi > 


backward branch of ^{E), W^{E) 


gi = 0, pi < 


forward branch of ^{E), ^^{E) 


= 0, gi > 


backward branch of ^{E), W^{E) 


= 0, gi < 


forward reactive spherical cylinder 
WAE) = W?(E) U W^(E) 


PiQi = 0, Pi,gi > 0, gi ^pi 


backward reactive spherical cylinder 
WtiE) = W,'iE)UWf^iE) 


Piqi = 0, pi,gi < 0, gi 7^ pi 


forward reaction path 


q2 = ■ ■ ■ = Qd = P2 = ■ ■ ■ = Pd = 0, Pi > 


backward reaction path 


q2 = ■ ■ ■ = Qd = P2 = ■ ■ ■ = Pd = 0, Pi < 



Table 1: Table of phase space surfaces influencing reaction dynamics and their representa- 
tions in normal form coordinates on an energy surface of energy greater than the energy 
of the saddle equilibrium point. 



However, all of these surfaces, and associated dynamical phenomena, are only 
"locally valid" in the neighbourhood U. The next step is to understand their influ- 
ence on the dynamics outside of U, i.e., their influence on the dynamics of reaction 
throughout phase space in the original coordinates (as opposed to the normal form 
coordinates). In order to do this we will need the normal form transformation con- 
structed in Section[2]and given in (|2.59p . to order N (where N is determined according 
to the desired accuracy following the discussion in Section 14. 5p . We rewrite p. 590 
below: 



A, — ^ ^U, 

(,r\ . . . , gir\pf\- ■ ■ = -(^^ = ° • • • ° ^'wM'')- (4-15) 

We refer to the original coordinates as the "physical coordinates" where reading 
from top to bottom, (j4.15|) describes the sequence of transformations from physical 
coordinates to normal form coordinates as follows. We translate the saddle-centre- 
• • • -centre equilibrium point to the origin, we "simplify" the linear part of Hamilton's 
equations, then we iteratively construct a sequence of nonlinear coordinate transfor- 
mations that successively "simplify" the order 3, 4, . . ., N terms of the Hamiltonian 
according to the algorithm described in Section [2j We can invert each of these trans- 
formations to return from the normal form coordinates to the physical coordinates. 

Computation of W^{E) and Wj{E): Our approach to computing the stable 
and unstable manifolds of a NHIM is, in principle, the same as for computing the 



stable and unstable manifolds of a hyperbolic trajectory (however, the practical im- 
plementation of the algorithm in higher dimensions is a different matter and one that 
deserves much more investigation). 

We describe the computation of Wj{E) as follows. 

• In the normal form coordinates, choose a distribution of initial conditions on 
the NHIM and displace these initial conditions "slightly" in the direction of the 
forward branch of W^{E) {pi = 0, qi = s > 0, e "small"). 

• Map these initial conditions back into the physical coordinates using the inverse 
of the normal form transformation. 

• Integrate the initial conditions forward in time using Hamilton's equations in 
the physical coordinates, for the desired length of time (typically determined 
by accuracy considerations) that will give the manifold of the desired "size". 
Since the initial conditions are in the unstable manifold they will leave the 
neighbourhood U in which the normal form transformation is valid (which is 
why we integrate them in the original coordinates with respect to the original 
equations of motion). 

The backward branch of W'^{E) can be computed in an analogous manner by dis- 
placing the initial conditions on the NHIM in the direction of the backward branch 
of WiE) ipi = 0, gi = e < 0, e "small"). 

Computation of W^{E) and Wj{E): The forward and backward branches of 
W'^{E) can be computed in an analogous fashion, except the initial conditions are 
integrated backward in time. 

Computation of the forward and backward reaction paths: Here the 
situation is, numerically, much simpler since we only have to integrate a trajectory. 
We consider the case of the forward reaction path. The backward reaction path 
is treated in the same way, after the obvious changes of sign for the appropriate 
quantities. 

Recalling that the dividing surface in normal form coordinates is given by qi = pi, 
the intersection of the forward reaction path with the dividing surface is given by 

q2 = ■ ■ ■ = qd = P2 = ■ ■ ■ = Pd = 0, 

ql =I,qi=pi> 0, with Kcnf(/, 0, . . . , 0) = ^ . (4.16) 

We transform this point in normal form coordinates into physical coordinates using 
the inverse of the transformations given in ()4.15p . Integrating this point forward in 
time using Hamilton's equations in the physical coordinates gives the forward reaction 
path immediately after passage through the dividing surface. Integrating the point 
backward in time gives the forward reaction path immediately before passage through 
the dividing surface. 

Computation of reactive volumes: Consider a region of the energy surface 
of some fixed energy E whose entrance and exit channels are associated with saddle- 
centre-- • • -centre equilibrium points. Near each such equilibrium point we can con- 
struct a dividing surface that a trajectory of energy E must cross in order to enter 
the region. Suppose that the region is compact and simply connected. An example 
is the phase space region associated with the potential well that corresponds to an 



isomer in an isomerization reaction |WBW04b] . It is then possible to give a formula 
for the energy surface volume corresponding to trajectories of the energy E that will 
leave that region of the energy surface. 

This formula is expressed in terms of the phase space flux across the dividing 
surfaces controlling access to this region of the energy surface and the corresponding 
"mean first passage times" of trajectories entering the region through the dividing 
surfaces. This theory is described in detail in [WBWOSal IWBWOSc] and here we 
just outline the results and show how the phase space structures discussed above in 
a region of the transition state are "globalized" to give this result. 

We consider an energy surface region to which entrance is only possible through 
a number of dividing surfaces, B'^'^~^^-{E) (i is the index for the number of forward 
dividing surfaces that control access to the region under consideration in the sense 
that trajectories initialized on this surface and integrated in forward time enter the 
region), and we compute the energy surface volume of reactive initial conditions, i.e., 
the initial conditions of trajectories that can leave the region under consideration 
through one of the dividing surfaces. The phase space transport theory described 
above is crucial for this computation as it allows us to define entrance and exit chan- 
nels uniquely in terms of dividing surfaces that have the property of "no recrossing 
of trajectories" and minimal directional flux. 

If the region under consideration is compact and connected it is a simple conse- 
quence of the Poincare recurrence theorem [Arn78j that reactive initial conditions in 
the region lie (up to a set of measure zero, or "zero volume") on trajectories which in 
the future escape from the region and in the past entered the region. Hence, for each 
point on a particular dividing surface hemisphere B'^^~^^-{E), there exists a time t 
(which depends on the point) for the trajectory starting at this point to spend in the 
region before it escapes through the same, or another, dividing surface. We define 
the mean passage time associated with B'^~^^{E) as. 




(4.17) 



Here we use the more concise language of diff'erential forms also used in Section [4.41 
to express the measure on the dividing surface over which we integrate the passage 
time. This measure is give hy Vt = uj'^~^ /{d — 1)!, where uj denotes the canonical 
symplectic two- form X]fc=i '^Pk ^ '^Qk- It then follows from arguments analogous to 
those that lead to the so called classical spectral theorem proven by Pollak in the 
context of bimolecular collisions |Pol81| . that the energy surface volume of reactive 
initial conditions in an energy surface region is given by 

Vreact(-E') = ^(i)enter;i(-E') fenter;i{E) , (4.18) 

i 

where the summation runs over all dividing surfaces -B^sli(^) controlling access to 
the region under consideration, and each entrance/exit channel contributes to the 
total reactive volume by the product of the associated mean passage time and the 
(directional) fiux, 

f enter AE) = [ ^. (4.19) 

The mean passage time for a given dividing surface hemisphere can be computed from 
a Monte Carlo sampling of that hemisphere. Performing such a sampling, uniformly 
with respect to the measure il, is straightforward in the normal form coordinates. 



The flux through a dividing surface hemisphere is also computed easily from the 
normal form as described in Sec. 14.41 The efficiency of this procedure has been 
demonstrated for concrete examples in |WBW05al IWBWOSc] 

Practical considerations: By their very definition, invariant manifolds consist 
of trajectories, and the common way of computing them, and visualizing them, that 
works well in low dimensions is to integrate a distribution of initial conditions lo- 
cated on the invariant manifold (hence, this illustrates the value of the normal form 
coordinates and transformation for locating appropriate initial conditions). In high 
dimensions there are numerical and algorithmic issues that have yet to be fully ad- 
dressed. How does one choose a mesh on a 2(i — 3 dimensional sphere? As this mesh 
evolves in time, how does one "refine" the mesh in such a way that the evolved mesh 
maintains the structure of the invariant manifold? 

4.7 The flux-flux autocorrelation function formalism for 
computing classical reaction probabilities 

In the chemistry literature (see |YT60l [M5T831 IMil98ap the accepted expression for 
the flux that goes in to the expression for the classical reaction rate is given by 

f{E)= [ [ S{E-H{q,p))F{q,p)P,{q,p) dqdp. (4.20) 

We want to explain the relation of this expression for the flux to the one derived 
in Section 14.41 We begin by explaining the dynamical significance of each function 
in (j4.20p . The function 5{E — H) restricts the integration to the energy surface of 
energy E under consideration. The remaining functions in the integral are defined 
on the basis of a dividing surface which is defined as the zero level set of a function 
s, i.e. the dividing surface is given by 

{{q,p)eR"' : s{q,p)=0}. (4.21) 

It is assumed that this surface divides the phase space into two components: a 
reactants component which has s{q,p) < and a products component which has 
siQ,p) > 0. In the chemistry literature s is usually a function of q only, i.e. "it is a 
dividing surface defined in configuration space." However, it is crucial to note that 
this restriction is not important. 

If we let Q denote the Heaviside function (which is zero if its argument is negative 
and one if its argument is positive) then the composition Q o s can be viewed as a 
characteristic function on phase space which vanishes on the reactants components 
and is identically one on the products component. The function F occurring in (|4.2U|) 
at a point {q,p) is then defined as the time derivative of 6os(<I>^(g,p)) at time t = 0, 
i.e., 

F{q,p) = ^^eos{'^'H{q,p))l^^ = 6{s{q,p)){s,H}{q,p), (4.22) 

where {•, •} again denotes the Poisson bracket. This means that F is a 6 function in 
s that is weighted by the scalar product between the gradient of the surface s and 
the Hamiltonian vector field Xh, 

F{q,p) = 6{s{q,p)){Vs{q,p),XH{q,p)) ■ (4.23) 



Due to the function 5{s) in F the integral (|4.2U|) is effectively restricted to the divid- 
ing surface (I4.2ip . or if we also take into acount the function 6{E — H), the integral 



()4.20p is effectively a (2d — 2) -dimensional integral over the intersection of the divid- 
ing surface (|4.21|) with the energy surface of energy E. It is not difficult to see that 
if we disregard the factor Pj- in (I4.20p . then the restriction of the resulting measure 
{s, H} dqdp to the intersection of the dividing surface with the energy surface agrees 
with the measure Q' that we defined in ()4.13p in Sec. 14.41 This implies that the ex- 
pression for the flux (j4.20p is invariant under symplectic coordinate transformations. 
The function Pj- in (j4.20p is defined as 

P,{q,p) = lim e(s(^*//(g,p)) , (4.24) 

t— +CXD 

which evaluates to one if the trajectory with initial conditions {q,p) has s{q{t),p{t)) > 
and hence proceeds to products for t — > oo and to zero otherwise. In this way 
the function Pj- in ()4.20p acts as a characteristic function on the intersection of the 
dividing surface with the energy surface. 
Equation (|4.2Up can be rewritten as 



f{E)= / CF{t)dt, (4.25) 
Jo 

where 

CF{t)=f [ 6{E-H{q,p))F{q,p)F{q{t),p{t)) dqdp, (4.26) 

which is referred to as the flux-flux autocorrelation function. This result is obtained 
from using the identity 

POO J 

Pr{q,p) = / -@os{'^'H{q,p))dt 

^?oo (4.27) 
F{^'H{q,p))dt, 



and changing the order of the time and phase space integrals. In (j4.27p it is tacitly 
assumed that Q{s{q,p)) = 0, which means that if we want to use the form of P^ given 
in (j4.27p in the integral (|4.20p then it is assumed that Q(s{q,p)) evaluates to zero 
on the dividing surface. This means that one assumes that a trajectory with initial 
condition at a point {q,p) on the dividing surface (j4.2ip still requires an infinitesimal 
time to actually cross the dividing surface (I4.2ip . i.e. more correctly (I4.25P should 
be 

/oo 
CF{t)dt. (4.28) 

We emphasised this point since it is important for understanding the time dependence 
of the function Cp to which we come back below. 

As stated in the chemistry literature (see e.g. |Mil98aj ) the equivalent expressions 
for the flux in (|4.20p and (|4.25p do not depend on the particular choice of the dividing 
surface. To see this recall that an arbitrarily chosen dividing surface will in general 
have the recrossing problem that we mentioned in Sec. 14.11 This means that there 
are either 

• "nonreactive recrossings" : nonreactive trajectories that cross the dividing sur- 
face, or 



"reactive recrossings" : reactive trajectories that cross the dividing surface more 
than once, 



or both. 

In fact reactive and nonreactive recrossings are independent, i.e. one can con- 
struct a dividing surface that only has nonreactive recrossings, or only has reactive 
recrossings or has both (or no recrossings at all like the dividing surface that we 
construct). Prom the definition of the function Pr in ()4.24p it is clear that those non- 
reactive recrossings that result from trajectories that approach the dividing surface 
from the side of reactants, cross the dividing surface ()4.2ip (two or an even number 
larger than two times) and return to the side of reactants do not contribute to the 
integral (I4.20p . In order to see that the factor Pj. in the expression for the flux in 
(j4.20p also takes care of nonreactive trajectories that approach the dividing surface 
from the products side and also of reactive recrossings one needs to note that F(q,p) 
takes into account in which direction a trajectory crosses the dividing surface: The 
sign of the scalar product between the Hamiltonian vector field and the gradient of 
the function s that defines the dividing surface depends on the direction in which 
the Hamiltonian vector field pierces the dividing surface (see (|4.23p ). In this way a 
family of nonreactive trajectories that approach the dividing surface from the prod- 
ucts side, crosses the dividing surface (j4.2ip (two or an even number larger than two 
times) and returns to the side of products will have a vanishing net contribution to 
the integral (j4.20p . Similarly, if a family of reactive trajectories crosses the dividing 
surface on its way from reactants to products n times (where n must be odd for the 
trajectories to be reactive) then the net contribution of the first n — 1 intersections 
of this family of trajectories to the integral (|4.20p is zero. This can be rigorously 
proven using the methods described in [W WQ4j but we omit the details here. 

The benefits that result from (|4.2Up formally not depending on the particular 
choice of the dividing surface are diminished by the fact that the implementation of 
the characteristic function Pj- is computationally very expensive. In practice (i.e., in 
numerical computations) one cannot carry out the integration of Hamilton's equa- 
tions to t = oo in order to evaluate according to (j4.24p . Instead one attempts to 
truncate the integration after a finite time Iq after which trajectories are assumed 
not to come back to the dividing surface. This is equivalent to assuming that the 
flux-flux autocorrelation function Ci?(t) is essentially zero for times t > to such that 
the integral in (j4.25p can be truncated at time to. A smaller time to required for 
this assumption to hold means that the amount of numerical computations required 
is reduced. This implies that some dividing surfaces are better suited for numerical 
computations than others |PM05j . but this is generally not known a priori. 

We note that our dividing surface is free of recrossings. In order to use expression 
(|4.2Up to get our result for the flux in (|4.1Up we deflne the function s according 
to s{q,p) = qi — Pi where {q,p) are the normal form coordinates that we used in 
Sec. 14.21 The delta function 6{E — H(q,p)) in the integral (|4.2ip then restricts the 
integration to the isoenergetic dividing surface that we constructed in Sec. 14.21 In our 
case Pj: simply needs to effectively restrict the integral (j4.2ip to the forward reactive 
hemisphere of our dividing surface. We therefore set 

Pr{q,p) = e{qi-pi). (4.29) 

In this way we recover the expression for the fiux that we have given in (|4.10p . It is 
crucial to note that in our case the evaluation of i-*r does not require the integration 
of Hamilton's equations and is therefore computationally much cheaper than using 
(j4.2Up with Pr defined according to (j4.24p for an arbitrarily chosen dividing surface. 
Equivalently, using the fact that in our case we have F = {H, P^} it is easy to see that 
the flux- fiux autocorrelation function Cpit) becomes the function 5{t) times our result 
for the flux given in (I4.10p . The time integration in ()4.28p (or in its corrected version 



()4.28p ) becomes trivial in our case. For an arbitrarily chosen dividing surface Cp 
will as a function of time gradually approach zero - in a monotonic or an oscillatory 
manner depending on the portions of reactive and nonreactive recrossings of the 
dividing surface (see e.g. [PM05j ). 



5 Quantum Reaction Dynamics and Cumula- 
tive Reaction Probabilities 

As described in the introduction, in this section we develop the quantum version of 
the classical reaction rate theory developed in Section HI We especially emphasize 
the roles of the classical and quantum normal forms. In particular, the classical 
coordinates in this section are the normal form coordinates. Moreover, we will see 
that the classical phase space structures that are realized through the classical normal 
form the "skeleton" on which the quantum dynamics evolves. 



5.1 Quantum normal form 

We consider a Hamilton operator whose principal symbol has an equilibrium point 
of saddle-centre-- •• -centre stability type. In Sec. 13.41 we have shown how such a 
Hamilton operator can be transformed to quantum normal form to any desired or- 
der N of its symbol by conjugating it with suitable unitary transformations. The 
resulting N^^ order quantum normal form H^^p is a polynomial of order in 
the operators 



.,^^d^ + 2j '^'= = -ydi + ^'^^' ^ = ^---d, (5.1) 



i.e., Hq^p is of the form 



-"QNF ~ -"■QNFV-' 5 -^2, • • • , 'Jd) 



Eq + \i + uj2j2 + ■ ■ ■ + ^d'Jd ch + higher order terms , 



(5.2) 



where c G B, is a constant and the higher order terms are of order greater than one 
and less than [A^/2] in the operators / and J^, k = 2, . . . ,d. 

From the structure of in (15.20 it follows that its eigenfunctions are prod- 

ucts of the eigenfunctions of the individual operators in (|5.1|) . This structure is the 
quantum manifestation of the integrability of the classical normal form described in 
Section 12.3.21 In the classical case integrability leads to a particular simple form of 
Hamilton's equations which provides a complete understanding of the phase space 
structure and dynamics in a neighborhood of the saddle-center-- - - -center equilibrium 
point. Similarly, we see that the quantum manifestation of classical integrability will 
lead to a simple structure for the corresponding quantum Hamilton operators in such 
a way that multidimensional problems are rendered "solvable". 

The operators are the Hamilton operators of one-dimensional harmonic oscil- 
lators (with unit frequency). Their eigenvalues are h{nk + 1/2), G INq, and the 
corresponding eigenfunctions are given by 

^-^^^^ = (vrn)V4V2».n.! ''"-(^)^"" ' ^'''^ 



where H^^. is the n}^ Hermite polynomial |AS65t ILLOlj . 

We will choose the eigenfunctions of I in such a way that their product with 
the harmonic oscillator eigenfunctions (j5.3|) give incoming and outgoing scattering 
wavefunctions of the system described by the Hamilton operator in (15. 2p . For clarity, 
we start with the one-dimensional case. 



5.2 Scattering states for one-dimensional systems 

The scattering states and S-matrix associated with a saddle equilibrium point in 
a one-dimensional system have been studied in |CP94al ICP94bl ICP99j and in the 
following we mainly follow their presentation. 

For one-dimensional systems a Hamilton operator in quantum normal form is a 
polynomial function of the operator / = —ih(^qd/dq + 1/2). The scattering states V'/ 
are the eigenfunctions of /, i.e., solutions of 

Hi{q) ^ -i^(g^ + ^) i^iiq) = I^iiq) (5-4) 



with eigenvalues I G R. Two solutions of this equation are given by 

V'/o;r(g) = e(-g)|grl/2+i^/\ 
V'/o;p(g) = e( q)\q\-'/'+''^\ 



(5.5) 



(5.^ 



where is the Heaviside function, and the index 'o' is for 'outgoing to' and 'r' and 
'p' are for 'reactants' and 'products', respectively. The motivation for this notation 
becomes clear from viewing the solutions (|5.5p as Lagrangian states, i.e., we rewrite 
them as 

^/o;r/p('Z) = ^7o;r/p('?)e''^^-/-('^)/^ , (5-6) 

where the amplitude and phase functions are given by 

Alo;r/piq) = ©(Tg)kr^/^ , ^/o;r/p(^) = I MqI , (5-7) 

respectively. This way we can associate the one-dimensional Lagrangian manifolds 

A/o;r = ^{q,P) = (^q, ^'^/o;r(g)^ ^ i^^'^^ : 9 < o| , 

A/ojp = I (?,?') = (^q, ^'^/o;p(g)^ ^ g) ' ^ ^ ^} 

with the states V'/o;r and tpio;p- From the presentation of A/o;r and A/o;p in Fig. [8] 
we see that for q — > — oo, ^/o;r is the outgoing state to reactants, and for q — > +00, 
ipio;p is the outgoing state to products. 

We define another set of eigenfunctions of / which will correspond to incoming 
states by requiring their momentum representations to be given by 

i^n-Ap) = V'/o;p(p) : '0/i;p(p) = ^/o;r(p) • (5-9) 

Here '*' denotes complex conjugation. The corresponding position representations 
are obtained from the Fourier transforms of ()5.9p giving 



(5.10) 



^PnAq) = [ ^n-M^^"" = I p-'/'-'^/^e^'^^' dp , 







i^n-Aq) = ^ / ^n-Av^^'" dp = ^ / (-p)-V2-i//^e^^P dp . 




Figure 8: Lagrangian manifolds A7o/i;i/p associated with the wavefunctions '0/o/i;r/p- The arrows indicate 
the Hamiltonian vector field generated by / = pq. 



The integrals in (jS.lOp are not absolutely convergent, but can be defined as oscillatory 
integrals. The motivation for defining incoming states according to Equation (|5.9p 
becomes clear from considering the stationary phase contributions to the integrals 
(jS.lOp . These come from the p satisfying 

^{-I\n\p\+qp)=^, (5.11) 
ap 

i.e., p = I /q, where p > for V/i;r and p < for V/i;?- This way we can associate 
with the incoming states the Lagrangian manifolds 



A/i;p = { iq,p) = [q,-] ■■ P<0} . 



(5.12) 



These manifolds are also shown in Fig. [8] and we see that for p — > +oo, ipn-r is an 
incoming state from reactants and for p — > — oo, ipn-p is an incoming state from 
products. 

In order to evaluate the integrals (jS.lOp we use the well known formula 



oo 



y^-^e-'^J' d2/ = e-"^°^r(z) . (5.13) 

This is valid for Re A; > 0, and we will use the analytic continuation to Re A; = 0, in 
which case the left hand side is defined as an oscillatory integral. We then obtain 

|l^e-i^l'^Viir(l-i£(-g)-i/2+i/A, q<0- ^ - ' 
This can be rewritten as 

V'/i;r = -^e-^^^°^r|^- - i-J {e^iipio-p - ie"^^V/o;r) ■ (5.15) 
In the same way we obtain 



For what follows in Sec. 15.71 it is useful to discuss how the eigenfunctions tpio-v/p 
and V'/iir/p s-re related to the more standard eigenfunctions of the operator / in the 
Q-representation that we introduced in Sec. 13.21 (see (j3.36p - (|3.38p ). 

The eigenvalue equation (15. 4|) then becomes 

hiiQ) = ( - f ^ - Iq') xdQ) = ixiiQ) . (5.17) 

Two solutions of this equation are given by 

where Di, is the parabolic cylinder function |AS65llLL01j . In fact, the eigenfunctions 
V'7i;r/p ^.re the images of Xi+/~ under the unitary transformation Ur that we defined 
in ()3.39p . or equivalently 

x/+ = ^;^/i;r, XI- = u: i^iup ■ (5-19) 

This relationship is discussed in great detail in |Chr03al IChr03b] where it is also 
shown that the pairs of eigenfunctions ipn-r/p^ '4'io;r/p and Xi+/- orthogonal and 
fulfill the completeness relations 



{rn./qW.Al') + rn-p{q)4^n;p{q')) dl = 6{q - q') , 

{rio-Aq)i^io-Aq') + rio,p{q)^io-M)) = s{q - q') , (5.20) 

/ {x%{Q)xi+{Q') + X'i-{Q)XI~{Q')) dl = 5{Q - Q') . 

5.3 S-matrix and transmission probability for one-dimensional 
systems 

The incoming and outgoing wavefunctions defined in Sec. 15.21 are not independent. 
Each solution ipj of (j5.4p can be written as a linear combination of V'/o;r/p i^Hiv/pi 

Ipl = ttplplo-p + arlplo;! , (5.21) 
Ipl = I3p1pn;p + MR;r ■ (5.22) 

These representations are connected by the S-matrix, 

We can read off the entries of the S-matrix from (I5.15P and (I5.16P and obtain 

Sm = ^e-'*'»«r(i - il) {-KT 'itx) . (5.24) 

Using the relation r(l/2 + iy)r(l/2 — \y) = 7r/cosh(7ry) it is easy to see that 
S{I)*S{I) = 1, i.e., S{I) is unitary. 

From the S-matrix we can determine the transmission coefficient 



and the reflection coefficient 



7^(/) = |5n(/)|^ = — J = . (5.26) 

e'^ft + e-'^ft 1 + e^'^fi 

As required we have T{I) + TZ{I) = 1. We see that the relevant scale is I /h. T tends 
to 1 if / > and to if K -h. 

We can generalise this now easily to operators -ffqNF = Kq^y{I)i where Kqnf is 
a polynomial function of /. In this case the incoming and outgoing states defined in 
Sec. 15.21 are also eigenfunctions of Hq^-p- We have 

(5.27) 

where E = Kq^-p{I) with / being the corresponding eigenvalue of /. The expression 
for the S-matrix in (j5.24p remains valid with / replaced by I{E) := Kq^p{E), where 
we have to assume that the energy is close enough to the equilibrium energy so 
that Kq^p{E) is invertible. We thus obtain the S-matrix for the scattering problem 
described by the Hamilton operator /?qnf = -f^QNF(-^) , 

S{E) = S{I{E)) . (5.28) 

The corresponding transmission coefficient is given by 

and similarly the reflection coefficient is given by R{E) = TZ{I{E)). This is a simple 
generalisation of the previous example. However, it is a very important result because 
we see that we can use the quantum normal form to compute the local S-matrix and 
the transmission and reflection coefficients to any desired order of the symbol of the 
Hamilton operator that describes the scattering problem. 



5.4 S-matrix and cumulative reaction probability for multi- 
dimensional systems 

We now consider the multi-dimensional case. In this case the Hamilton operator in 
quantum normal form is given by -f^QNF = Kq^p{I, J2, ■ ■ ■ , Jd), where -fCgNF is a 
polynomial function, and = {—h'^dg^ + (?|)/2, k = 2,. . . ,d, are one-dimensional 
harmonic oscillators. Let V'nfc) '^fc £ I^O) be the n^*^ harmonic oscillator eigenfunction 
(I53D, i.e., 

JfcVn, = Kuk + l/2)V'n, . (5.30) 

Then the incoming and outgoing scattering states are given by 

(5.31) 

,naca) o;r/p 

where Ugca = (n2, . . . , nd) G INg Ms a (d - l)-dimensional vector of scattering quan- 
tum numbers. 

The S-matrix connecting incoming to outgoing states is then block-diagonal with 



(5.32) 



where (5nsca,msca is the multi-dimensional Kronecker symbol, S(I) is given by (I5.24p 
and /nsca(^) is determined by 



i^QNF(/n_(i^), Kn2 + 1/2), . . . , h{nd + 1/2)) = E 



(5.33) 



We will assume that this equation has a unique solution In^^^ (E) , which is guaranteed 
if the energy is close enough to the equilibrium energy since i^QNF starts linearly in 
the actions, see (15. 2p . 

We can now define the transition matrix T as the diagonal sub-block of the S- 
matrix which has the (1, 2)-components of the matrices in (I5.32p on the diagonal, 
i.e., 



sea 1 ' ' I'sca 



cai " I'sca 



5l,2(/n_(i^)) 



1 + exp 



27r 



h 



-1 



(5.34) 



The cumulative reaction probability N{E) is then defined as (see, e.g., |Mil98a| ) 

N{E) = Tr T{E)T{E)'' . (5.35) 

Using (j5.34p we thus get 



N{E)=Y,Tr. 



'sca ) '"-sea 



E 



1 + exp 



27r- 



h 



(5.36) 



The cumulative reaction probability N{E) is the quantum analogue of the classical 
flux f{E) or, more precisely, of the dimensionless quantity N-y\ieyi{E) = f{E) / 
that we defined in Equation (14. lip in Sec. 14.41 To see this let us consider N{E) in 
the semiclassical limit /i — > 0. To this end first note that 



1 + exp 



27rl/hj 



-1 



e(/) ash^o, 



(5.37) 



where is the Heaviside function. This means that the transmission coefficients 
^rasca,nsca(-^) i'^ (|5.36p arc essentially characteristic functions, i.e., in the semiclassical 
limit, r„_,„_(£;) is or 1 if the solution of K(/„_, ;i(n2 + l/2), . . . , h{nd+l/2)) = E 
for /nsca is negative or positive, respectively. This way the cumulative reaction proba- 
bility can be considered to be a counting function. For a given energy E, it counts how 
many of the solutions of the equations i^QNF (-^nsca i ^(^^2+1/2), . . . , h{n(i+l/2)) = 
E with scattering quantum numbers risen = {^2, • • • , nd) G Mq"^ are positive: 



NiE) ^ > 



KQNF{In,,,,Hn2 + ^),..., h{nd + ^)) 



E, n,ca G , 



(5.38) 



as — > 0. In other words, N{E) can be considered to count the number of open 
'transmission channels', where a transmission channel with quantum numbers risca is 
open if the corresponding transmission coefficient 2^rasca,nsca(-^) is close to 1. 

We can interpret N(E) graphically as the number of grid points (/i(n2+l/2), . . . , h{n 
1/2)) in the space of {J2,---,Jd) S [0,oo)'^~^ that are enclosed by the contour 
-f^QNF(0, J2, . . . , Jd) = E, see Fig. O The number of grid points is approximately 
given by the volume in the space of ( J2, . . . , Jd) £ [0, oo)'^"^ enclosed by -K'qnf(0, J2, • • • , 
E divided by h'^~^. Using the fact that for — > 0, -fiTqNF becomes the function -fCcNF 
which gives the classical energy as a function of the classical integrals (/, J2, . . . , Jd) 




Fi gure 9: (a) Lines {I,h(n2 + 1/2), .. . ^h{nd + 1/2)), / € H, £ INg, k — 2, . . . , d, in the space 
(/, J2, . . . , Jd) G R X [0,oo)''~^ for d — 3 and their intersections with the surface ifQNF(/,j2,J3) — ^- (b) 
Grid points {h{n2 + 1/2), . . . , h{nd + 1/2)) in the space (J2, . . . , Jd) for d — i. The blue line marks the 
contour iCQNF(0, J2, ■ • • , Jrf) = i?- In this plot only the scattering states for which the quantum numbers 
(712,77.3) have the values (0,0), (0, 1), (1,0) or (1, 1) correspond to "open transmission channels", see text. 



we find that the volume in the space of (J2, . . . , Jd) enclosed by ETcnfCO, J2, • • • , Jd) = 
E is given by the classical flux f{E) divided by (27r)'^~^, see (|4.10p in Sec. l4.4[ and the 
cumulative reaction probability N{E) is thus approximately given by N-sj^cyiiE) = 
f{E)/(27rh)'^~^ defined in (j4.11|) in Sec. 14. 4[ This way we verified our statement in 
Sec. l4.4l that A'weyi(-E') gives the mean number of open transmission channels. In fact, 
as mentioned in Sec. 14.41 the classical flux f{E) can be considered to be the phase 
space volume enclosed by the energy contour of energy E of the invariant subsystem 
which has one degree of freedom less than the full scattering system and which as 
the so called activated complex is located between reactants and products. A'weyi(-E') 
counts how many elementary quantum cells of volume {27rh)'^~^ fit into this phase 
space volume and this way gives the Weyl approximation of the cumulative reaction 
probability N{E). 

It is important to note here that like the flux in the classical case the cumula- 
tive reaction probability is determined by local properties of the Hamilton operator 
embodied in its symbol in the neighbourhood of the equilibrium point only. All one 
needs to know is the quantum normal form, which enters through the relation (15.330 
and which determines In^cS^)- 

5.5 Distribution of the scattering states in phase space 

At the end of the previous section we have seen how the cumulative reaction prob- 
ability is related to the classical flux. In this section we want to further investigate 
the quantum classical correspondence by studying the distribution of the scattering 
states in phase space and relating these distributions to the classical phase space 
structures that control classical reaction dynamics as discussed in Sections 14.11 and 

The standard tool to describe the phase space distribution of a wavefunction is the 
Wigner function, but since the scattering wavefunctions are not square integrable the 
Wigner functions will be distributions. Therefore it is more convenient to study the 
phase space distribution in terms of their Husimi representation which is obtained 
from projecting the scattering states onto a coherent state basis (see |Har881 FBalQS] ) 



and this way leads to smooth functions. For a point {qo,Po) ^ >^ we define a 
coherent state with wavefunction 



V',o,Po(9) = ^-^ei«PO'^)-<«0'Po)/2)e-2k('?-5o,.-,o> . (5.39) 

This wavefunction is concentrated around q = Qo and its Fourier transform, i.e. 
its momentum representation, is concentrated around p = Po- In phase space the 
coherent state (|5.39p is thus concentrated around (qo^Po)- The Husimi function of a 
state '0 is now defined by the modulus square of the projection onto a coherent state, 

H^iQ,P)-=J^\{i^P,,,^)\'- (5.40) 

It has the important property that the expectation value of an operator Op [A] with 
respect to a state is given by 



{ij,Op[A]^) = / / A{q,p)H^{q,p) dqdp+Oih) . (5.41) 

Furthermore, we have H^(q,p) > 0, i.e., the Husimi function can be considered to be 
a probability density on phase space and describes how a quantum state is distributed 
in phase space. 

The Husimi functions of the scattering states V'(/,nsca) i/o;r/p inherit the product 
structure (|5.3ip . i.e. we have 

'^V'(/,„eca)i/o:r/p(9i' ■ ■ ■ , Qd , Pi , ■ ■ ■ , Pd) = H^j,^^,^/^iqi,Pi)H^„^ (^2,^2) " " " H^„^ {qd,Pd) . 

(5.42) 

The Husimi functions of the eigenfunctions tpnk of the one-dimensional harmonic 
oscillators Jk are well known (see, e.g.. [KMW97j ). 

The first three of these Husimi functions are shown in Fig. [TOj They are concentrated 
on the circles pf. + q^. = 2nkh and have an n^-fold zero at the origin. 

The computation of the Husimi functions for the one-dimensional scattering states 
4'io;v/-p in (15. 5p can be found in |NV97] where it is shown that for the linear combi- 
nation 

V'"'^ = aV/o;p + /3V'/o;r , a, /3 G C , (5.44) 



one gets 



^ l-Kh cosn.{TT 1 /h) 



(5.45) 



where Dy again denotes the parabolic cylinder function |AS65j . Fig. [TOl shows con- 
tour plots of the Husimi representation of the state tpri;r for different values of the 
eigenvalue I. Here a and /? in (j5.44p are determined from (j5.15p . In accordance with 
the classical dynamics where trajectories with I < are non-reactive and trajectories 
with / > are reactive, most of the state '(/'/i;r is refiected to the reactants side for 
/ < while it is mainly transmitted to the products side for / > 0. The border- 
line case between these two situations is given by I = 0. Here the state is localised 




■:<.-?■ Ij 14- 



Figure 10: Contour plots of the harmonic oscillator Husimi functions H.^^^ in the ((7/i;,p/i;)-plane for 
rife = (a) Tik — 1 (b) and = 2 (c), and contour plots of the Husimi fimctions H^^^,^ in the (gi,pi)-plane 
for / = — 1 (d) / = (e) and / = 1 (f). Red corresponds to low values; blue corresponds to high values. 
In (a)-(c) the spacing between the values of the contourlines is decreasing exponentially, {h = 0.1 .) 



in phase space at the hyperbolic equilibrium point with ridges along the reactants 
branches of the stable and unstable manifold and the products branch of the unstable 
manifold. 

Fig. [10] indicates that the Husimi functions of the scattering states ipn-r are lo- 
calised on the Lagrangian manifolds 

A-(/,nsca)i/o;r/p = A/i/o;r/p X ^n2 X • • • X , (5.46) 

where the A/i/o.i./p are defined in ()5.8p and (15.120 . and 

An, = {{qk,Pk)&^^ ■■ ql+pl = 2hnk}, k = 2,...,d, (5.47) 

are the Lagrangian manifolds associated with one-dimensional harmonic oscillator 
eigenf unctions. Quantum mechanics thus picks out those Lagrangian manifolds 
^fj^ Jd foliating the classical phase space (see Sec. 14. 3p for which the actions, 
J2, . . . , Jd, fulfill Bohr-Sommerfeld quantisation conditions. More precisely we find 
that the outgoing scattering states ipi-o;r/p localised on the the Lagrangian man- 
ifolds 

^(/,nsca)o;r = ^1 ,hn2,...,hnd ' ^g-j 
^(/.risca) o;p = ^I,hn2,...,hnci ' 

and the incoming scattering states 'i/'/;i;r/p localised on the Lagrangian manifolds 

A _ / ^IM2,-Md ' ^ >^ 

^»-(/,nsca)i;p 1 A+ T ^ Ci ' 

I ^^I,hn2,-,hnd ' ^ 

The projection of the Lagrangian manifolds A{i^ns„s.)i/o\r/p to the centre planes 
ilkiPk), k = 2, ... ,d, is thus restricted to the discrete circles p1+ql = 2nkh, Uk E INq. 




Figure 11: Projections of the Lagrangian manifolds A(/ „^^^)i.i. defined in Equation (|5.46p to the normal 
form coordinate planes for the same setup as in Fig.[9l The scattering quantum numbers are ngca — {'^2, n^) 
with < 77.2,713 < 3. For the values (0,0), (0,1), (1,0) and (1,1) of the quantum numbers (71,2,77^3), the 
Lagrangian manifolds A^in ^ are contained in the energy surface volume (green region) enclosed by the 

forward reactive spherical cylinder Wf{E) defined in Sec. HI For the other values of the quantum numbers 
the Lagrangian manifolds A(7 „^^^) i-r are located in the reactants component of the energy surface. 



If we fix the total energy E then this also entails a discretisation of the projection of 
the manifolds (|5.46|) to the saddle plane {qi,pi) since the eigenvalue / needs to satisfy 
the energy equation Kq^^p^I, h{n2 + 1/2), . . . , fi{nd + 1/2)) = E. For the Lagrangian 
manifold ^(^i^nsca)i;r ^^^^ depicted in Fig. [TTl Depending on whether / is positive 
or negative the Lagrangian manifold A(/ „^^_^)i.i. is either located inside or outside of 
the energy surface volume enclosed by the forward reactive spherical cylinder Wj {E) 
defined in Sec. HI and hence is either composed of reactive or nonreactive trajectories 
of the classical dynamics. From our discussion at the end of Sec. 15.41 it then follows 
that the cumulative reaction probability N{E) is approximately given by the total 
number of Lagrangian manifolds A(/ y,. which, for scattering quantum numbers 
"-sea = in2, • • • , rid) £ lNg~^, are located inside of the energy surface volume enclosed 
by Wf{E). 

5.6 The global S-matrix 

It is important to emphasize again that, so far, our approach to quantum reaction 
dynamics has been local, i.e., it is derived completely from the properties of the 
quantum normal form that is valid in the neighborhood of the saddle-centre-- •• - 
centre equilibrium point. The property of the resulting S-matrix in (j5.32p being 
block-diagonal reflects the fact that the quantum normal form is integrable in the 
sense that the basis of scattering states can be chosen in the product form (j5.3ip . In 
a different basis the matrix will lose this feature, and phenomena like mode mixing 
are related to how other incoming and outgoing scattering states are related to this 
special basis. It is natural to embed the study of this phenomenon in a study of the 
global dynamics which we will describe in this section. The global formalism is in 
particular required in order to compute general state-to-state reaction rates. 

Let us start by describing the scattering or reaction process in classical mechanics 
by using Poincare sections. Recall that a Poincare section at energy E is given by a 
smooth hypersurface T,{E) of the energy surface with energy E which is transversal to 
the flow (S(ii^) is allowed to have several components). If we have two such Poincare 
sections Tii{E) and 'E2{E) such that all the flow lines intersecting T,i(E) intersect at 
a later time 'E2{E), too, then moving along the flow from T,i{E) to S2(i?) defines a 
Poincare map 

:5]i(E) ^S2(^) . (5.50) 
Such Poincare maps can be composed. If T,3(E) is another Poincare section which lies 



behind 5]2(i?) in the sense that the flow lines that intersect Y^2iE) also intersect '^^{E) 
at a later time, and if P^^''^\E) : T,2{E) Tj^{E) is the corresponding Poincare map, 
then the Poincare map 

p(3,i)(i^):Si(i^)^S3(i^) (5.51) 

is given by 

p(3,l) = p{3,2) ^ p(2,l) _ (5,52) 

Using this construction we can describe transport through phase space regions 
by a sequence of maps. Given some Poincare section Sinitiai(-£') located in the area 
of initial points in the reactants region where we prepare the system and a Poincare 
section Efinai(^) in the products region where we measure the outcome, a succession 
of Poincare maps 

Sinitial(^) ^ Si(^) ^ ^2{E) ^ > ^RnAE) (5.53) 

tells us how the initial points are transported through the systemic 

The advantage of subdividing the flow into a sequence of maps lies in the fact 
that different regions in phase space might need different techniques to compute the 
flow. In our case of interest Poincare sections can be constructed to the products and 
reactants side of a saddle-centre-- • • -centre equilibrium point. The dynamics 'across' 
this equilibrium point can then be described by the normal form while the dynamics 
between neighbourhoods of different saddle points can be obtained from integrating 
the original equations of motions |Cre041 ICreOSi I WB W05b] . Moreover, the phase 
space structures obtained from the local normal form can be "globalized" following 
the discussion in Section 14.61 

A similar procedure can be developed in the quantum case. The Poincare maps 

P^j'^HE) : J:,{E) ^ J^jiE) (5.54) 

are symplectic maps, and as such can be quantised using the theory of Fourier integral 
operators. The quantisations will be unitary operators which we interpret as local 
S-matrices, 

5(^'^)(i?):L|^(^)-L|^,(^) , (5.55) 

where L'^(^^^ is a Hilbert space obtained by geometric quantisation of ^{E), see, e.g., 
[KirOlj . This is similar to the quantisation developed in Bog92 . As in classical 



dynamics we can compose these matrices to obtain a global S-matrix 

^(final,initial) _ ^{final,n) (^)^(n,n-l) _ _ _ ^(l.initial) (5.56) 

which tells us how initial states in Li. are transformed into final states in 

^Sfl i{E) • '^^^ reasons for introducing this splitting of the S-matrix are the same as in 
the classical case. We can employ different techniques for computing the S-matrices 
according to different local properties of the system. Near equilibrium points the 
dynamics can be described by the quantum normal form we developed in this paper. 
Notice that the neighbourhoods of the saddle-centre- • • -centre equilibriuml points 
are the regions where we expect quantum effects to be of most importance due 
to partial reflection at and tunnelling through the barriers associated with saddle 
points. The quantum transport between neighbourhoods of different equilibrium 
points can be described by a standard van Vleck type formalisms, using, e.g, initial 
value representations (IVRs) which are very common in theoretical chemistry (see, 
e.g., |Mil98at iMUOSb] for references). 



We here ignore the difficulties involved in constructing global Poincare sections (see, e.g., |DW95| ): 
we assume that the sequence of Poincare sections (|5.53p is intersected transversally by the trajectories 
with initial points from a suitable open subset in the reactants region. 



5.7 The flux-flux autocorrelation function formalism to 
compute quantum reaction probabilities 



The main approach to compute quantum mechanical reaction rates that is most 
heavily pursued in the chemistry literature is the quantum version of the flux-flux 
autocorrelation function formalism that we reviewed in Sec I4.7[ This approach was 
developed by Miller and others (see [YTOni IMST831 IMil98aj ^ and in the following 
we will mainly follow their presentation. We will see that the cumulative reaction 
probability N(E) is the quantum mechanical flux through a dividing surface and 
hence is the analogue of the classical flux. The goals of this section are twofold. 
Firstly, we will show that we recover our result for the cumulative reaction probability 
in (j5.36p when we evaluate the quantum flux-flux autocorrelation function expression 
for the cumulative reaction probability N{E) in terms of the quantum normal form 
and for our choice of the dividing surface that we discussed in Sec. HI This way will see 
that the flux-flux autocorrelation function formalism and our result for the cumulative 
reaction probability are formally equivalent and hence, our result for N{E) can be 
viewed as a quantum mechanical flux through a dividing surface. Secondly, we will 
argue that, like in the classical case, the application of the flux- flux autocorrelation 
formalism in its original form, which does not depend on the specific choice of a 
dividing surface, is computationally much more expensive than our quantum normal 
form approach. 

Following |YT60[ [MST831 IMil98aj . a quantisation of the flux-flux autocorrelation 
function formalism in Sec. 14.71 or more precisely of the dimensionless quantity 

NweAE) = f{E)/{27Thf-' = 27Th [ [ 5{E - H)FP, (5.57) 

is obtained by replacing the classical phase space integral in (I5.57P by the trace of 
the associated operators in the form 

N{E) = 2TThTv6{E - H)FPr. (5.58) 

Following the quantum classical correspondence principle the operator F is obtained 
from its classical counterpart F by replacing the Poisson bracket in the classical 
expression F = {Q{s),H} by the corresponding commutator to give 

F = -^[e(^),H]. (5.59) 

Here Q{s) is a quantisation (to which we will come back below) of the composition of 
the Heaviside function with a function s that defines the dividing surface according 
to s{q,p) = as discussed in Sec. 14.71 Similarly, the quantisation of the projection 
function = limt^oo 6(s(^*)) in (j4.24|) is given by the operator 

= lim e^^*G(s)e-^^* . (5.60) 

The application of P^ to a state tp is thus obtained from taking the limit f ^ cxd in 
the process of letting the time evolution operator, exp(— -^i^t), act on -0 for the time 

t, then apply 0(s) to determine whether ip has evolved to products after time t (see 
below for the details), and then evolve the state ip backward in time by applying the 
inverse of the time evolution operator, exp(^//t). In fact, the operator P^ is given 

by the limit t ^ co of the Heisenberg picture of the operator G(s). 



Using 



P,= I A (ei^*e(5e-x^*) dt 







oo 



(5.61) 



we find that analogously to (j4.25p the cumulative reaction probability can be rewrit- 
ten as an autocorrelation functior^l: 

POO 

N{E) = 27rh Cp{t)dt, (5.62) 
Jo 

where 

Cp{t) = TV (5(S - H)FeT^"^Fe-^"^ . (5.63) 

We illustrate the application of the flux-flux autocorrelation function formalism 
in the following sections. 

5.7.1 Example: ID parabolic barrier 

As a first example we consider a one-dimensional system and a surface defined ac- 
cording to s{q,p) = q — qo = 0. In the position representation the quantisation of the 
function Q{s) is then defined by its action on a wavefunction ip{q) according to 

&is}^iq) = Q{q - qoMq) • (5.64) 

A state thus is an eigenfunction with eigenvalue 1 of the operator Pj- if its wave- 
function V(^) is concentrated in g > go if evolved forward in time to time t = oo. 
Likewise, ip is an eigenfunction with eigenvalue of the operator Pj. if its wavefunc- 
tion ip{q) is concentrated in g < go if evolved forward in time to time t = oo. For 
a Hamilton operator of type 'kinetic plus potential', H = + V{q), with the 

quantisation of the operators g and p given in (|3.38|) . the operator F becomes 



(5.65) 



2m 

For the expectation value of F with respect to a state we thus ge10 

= -il-{^l^*{qo)ij'{qo) - i^'*{qoMqo)) , (5.66) 

where the primes denote the derivatives. This agrees with the standard definition of 
the quantum probability current density that can be found in any quantum mechanics 
textbook (see, e.g., |LL01] ). 



^•^Formally Eg. 15.611 still contains a term Q{s). But this term will give no contribution to N{E) for the 
same reason as in the classical flux- flux autocorrelation formalism (see the discussion after ()4.27p '). In the 

examples below this can be seen explicitly since we define the operator 0{s) in normal form coordinates 
as a multiplication operator by a characteristic function. Then the same reasoning as in the classical case 
applies. 

^^In the following it will be notationally more convenient to use the Dirac notation for scalar products. 
Here is the same as ('0, Atp) for any operator A and state ip. 



To make the example more concrete we consider a parabolic barrier described by 
the Hamilton operator 

The spectrum of H is 11. We choose energy eigenfunctions ipE± such that they 
correspond to wavefunctions moving in positive and negative q direction, respectively, 
i.e., besides 

H^E±=E^E± (5.68) 

we have 

PrlljE + = V'E + , Pri^E - = . (5.69) 

For the trace (|5.58|) to be well defined we need to require that the states ipE± are 
normalised in such a way that they satisfy the completeness relation 



{rE + {Q)^E + {q) + rE-iQ)i^E-{q')) dE = 6{q - q') . (5.70) 
The eigenfunctions TpE± having the properities ()5.69p and (|5.70p are given by 



1 / m \V4 fl E\ / /2mA \ 

where again denotes the parabolic cylinder function |AS65j . In fact, the wave- 
functions ipE± can be obtained from a suitable scaling of the wavefunctions xi± that 
we defined in (j5.18p and which satisfy the completeness relations ()5.20p . For ipE±, 
we have 

- -^l^irE±i^'E± - i^'E±^E±) = ±^ i^^-L/(A.) ' (5-72) 

and hence using ()5.66p and ()5.69p we get for the cumulative reaction probability, 

N(E) = 27r/iTr 5{E - H)FP, 

= 27rh [ {{^JE' + \S{E - H)FP,\i;E' +) + i^E' -\SiE - H)FP,\^E' -)) dE' 
Ju 

= 27rh [ 5[E - E'){ijE' + \F\ilJE'+)dE' 



1 



1 _)_ Q-2nE/{\h) ' 

(5.73) 

which is the exact quantum mechanical reflection coefficient for a parabolic barrier 
[LLOT] . 

We now want to repeat the calculation above by inserting for H the quantum 
normal form of the parabolic barrier in (j5.58p . This will show two things. Firstly, 
this will lead to our result for the cumulative reaction probability N{E) that we have 
given in (I5.36P (which for the one-dimensional case reduces the reflection coefficient 
derived in Sec. l5.3p . Secondly, we will see that our result agrees with N{E) in (j5.73p . 
i.e., our result for N{E) in terms of the quantum normal form is exact for parabolic 
barriers. 

From our discussion in Sec. 13. Sl it follows that the quantum normal form of (I5.67P 
is given by 

^QNF = i^QNF(/) = A/. (5.74) 



In order to evaluate (I5.58|) for our dividing surface which in terms of the normal form 
coordinates is given by s{q,p) = q — p = (see Sec. 14. 2p it is convenient to work with 
the rotated coordinates 

{Q,P) = ^iq-p,q + p). (5.75) 
The Q representation of the operator 0(s) is then defined analogously to (|5.64p . i.e., 

G(5v(Q) = QiQMQ) ■ (5.76) 

As we have seen in the example of the application of Lemma [H] (exact Egorov) in 
Sec. 13.21 the Q representation of the operator I reads 

(see Equation ()3.38p ). In Sec. 15.21 we showed that the eigenfunction of ()5.77p are 
given by xi± defined in (I5.18p . In fact, the eigenfunctions xi± formally agree with 
the eigenfunctions in (|5.7ip if m and A are replaced by 1, and E is replaced by 
/. Analogously to (I5.65|) we have 

- ^[©S, I] = limQ) + m)p) , (5.78) 

and for an arbitrary state ip, 

(V'l - ^m^),m = -i|(v*(o)V''(o) - i^'*mm) ■ (5.79) 

Evaluating this expression for the eigenfunctions xi± we get 

- i-{x*i±x'i± - X7±XI±) = ±^ ^ ^ • (5.80) 

Using this result and the fact that xi+ and xi- moving in positive and negative 
Q direction and hence are eigenfunction of with eigenvalues 1 and 0, respectively, 
we get 

NiE) = 2TThTr 5{E - Kc^^Y{i))FP, 

= 27rh [ i{xi + \S{E - KciMl))FPr\xi+) + {xi-\S{E - KQMl))FPr\xi -)) dl 

= 27Th [ 6{E-KQ^F{I))X{xi + \-UQ{^)J]\xi+)dI 
1 



1 _)_ Q-2TrE/{Xh) ■ 

(5.81) 

This formally agrees with the expression for N{E) that we have given in ()5.36p and 
also with the exact result in (j5.73p . i.e., our quantum normal form computation of 
N(E) is exact for parabolic barriers. 



5.7.2 Example: General barriers in ID 

Let us now use the quantum normal form in the flux-flux autocorrelation formalism in 
the more general case of a one-dimensional system with a Hamilton operator whose 
principal symbol has a saddle equilibrium point but is not necessarily quadratic. 



Like in the previous section we again work in the Q representation, i.e., our dividing 
surface is defined by s{Q^ P) = Q = 0, and the operators Q{s) and / are defined by 
()5.76p and ()5.77p . respectively. In order to evaluate (I5.58|) for a general Hamilton 
operator in quantum normal form, -ffqNF = Kq^p{I), where /Cqnf(-^) is a polynomial 
in /, we use that for n G IN, we have 

rt-l 

[e5),/"] = J2r''-''-'m^)J]i''- (5.82) 

A:=0 

This can be shown by direct calculation. For the eigenfunction xi± of I we thus 
have 

{Xl±\ms),h\xi±) = {xl±\ms),i]\xi±)nr-\ (5.83) 

and hence 

{xi±\m),KQMmXi±) = {Xi±\m^)j]\xi±)^^^f^. (5.84) 

Using this together with (|5.79|) and (j5.8U|) we find that the cumulative reaction prob- 
ability is given by 

N{E) = 27Th iixi + m - KQMi))FP,\xi+) + ixi^m - AQNF(/))FPr|x/-)) 

= 2nh jj{E- i^QNF {I)){xi + \-'jms)J]\xi+) ^^^^ dl 
1 



1 + Q-27Tl{E)/h ' 

(5.85) 

where I{E) is the solution of E = KQf^p{I{E)), and we have assumed that there is 
only one such solution (compare with the remark after (|5.33p ). We thus recover our 
result for N(E) that we have given in ()5.36p . 

5.7.3 Example: General barriers in arbitrary dimensions 

We now consider the d-dimensional case with a Hamilton operator in quantum nor- 
mal form given by -f^QNF = -^qnf(-^5 J2, ■ ■ ■ , Jd)- Again we work in the Q representa- 
tion in terms of which our dividing surface is defined as s{Qi, . . . , Qd, Pi, ... , Pd) = 
Qi = 0. The quantisation of 0(s) is then defined by its action on a wavefunction 
ip{Qi, ■ ■ ■ , Qd) according to 

ei^mQi, ■■■,Qd) = e(Qi)^(Qi, ...,Qd). (5.86) 
The Q representation of the incoming eigenfunctions ()5.3ip is given by 



X(/,neca)i;r(Ql' ■■■,Qd) ■■= XI + {Ql)'>Pn2{Q2) ' ' ' i^n^Qd) , 
X{7,n,ca)i;p(Ql' ■■■,Qd) ■■= Xl^{Ql)'4'n.2{Q2) ' ' ' i'nAQd) 



(5.87) 



with I G K, and scattering quantum numbers ngca = (^^2, . . . ,nd) G Mq ^ . It then 
follows from the one-dimensional case discussed in the previous section that 

(X(/,neca)i;rl[0(s)>^QNF(i', ^2, . . . , ^d)] I X(/,n,ca) i;r) 

^ d 1 1 (5-88) 

= (X(/,n,ca)i;rl [0(s), -^]|X{/,nsca)i;r) qjKq^f{I, h{n2 + -),..., h{nd + -)) 



(see Equation (I5.84p ). Using the completeness of the states X(7',nsca) i;r/p '^^ find for 
the cumulative reaction probability, 



N{E) = 27rh Yl [ ((X(/,n_)i;r|5(^ " i^QNpU, J2, • • • , Jd))i^^^r|X(7,n...) i;r) 



+ {XiI,n,,,)i;p\S{E - Kqnf(/, J2, • • • , Jd))FPr\X{I,n,,,)i;p)) d/ 



27rh 



J2 [ S{E- i^QNF(/, ^("2 + • • • > + ^))) X 



(X(/,n,ca)i;rl - T[Q{s),KQT^F{i, J2, ■ ■ ■ , ^d)] |X(/,n,ca) i;r» d/ 



1+ exp — 27r 



h 
h 



-1 



(5.89) 



where In^^A^) solves i^QNF(/, ^(n2+l/2), . . . , h{nd+l/2)) = E for 



INq ^, and we assume there is only one such solution (compare, again, with the re- 
mark after ()5.33p ). We thus recover our result in ()5.36p . 

Though we showed that if the flux-flux autocorrelation function formalism is 
evaluated in terms of the quantum normal form then it reproduces our results for 
the cumulative reaction probability that we developed in Sec. 15.41 it is important to 
point out the computational differences between the flux-flux autocorrelation function 
formalism in its original form and the quantum normal form approach to compute 
cumulative reaction probabilities. The main problem with the implementation of 
the flux-flux autocorrelation function formalism is the occurrence of the projection 
operator in the trace in (j5.58p . The presence of the operator is crucial in order 
to ensure that only states that evolve from reactants to products contribute to the 
trace in (j5.58p . The extraction of this information for an arbitrarily chosen dividing 
surface and without any insight into the quantum dynamics requires one to look at 
the full time evolution of states as embodied in the definition of the operator P^ in 
(j5.60p . Though various techniques like Monte Carlo path integration and initial value 
representation (IVR) [Mil98a, Mil98b] have been developed in order to solve this time 
evolution problem that is involved in the evaluation of the trace in (|5.58p due to the 
presence of P^ it remains a formidable numerical task to apply (I5.58P to specific 
systems. In contrast to this, the computation of the cumulative reaction probability 
from the quantum normal form does not involve the solution of a time evolution 
problem. The reason for this is that the quantum normal form yields an unfolding 
of the quantum dynamics in the reaction region. As a result the S-matrix expressed 
in terms of the corresponding scattering states is diagonal, i.e., the scattering states 
can be immediately classified and the reaction probabilities can be immediatedly 
determined without explicitly looking at the time evolution. The numerical effort 
to implement and evaluate the quantum normal form is comparable to the classical 
normal form computation described in Sections [2] and [H In Sec. [7] we will illustrate 
the efficiency of the quantum normal form computation of the cumulative reaction 
probability for several concrete examples. 



(^2, 



6 Quantum Resonances 



In this section we consider quantum resonances and the corresponding resonance 
states. The role of quantum resonances in the context of chemical reactions has 
been studied for the first time explicitly in the chemistry literature by Friedman and 
Truhlar [FT91] and Miller |SM91] . The Quantum resonances are viewed as another 
imprint of the activated complex in addition to the quantisation of the cumulative 
reaction probability discussed in the previous section, Sec. O Recent developments in 
high resolution spectroscopic techniques allow one to probe the dynmaics of quantum 
mechanical reactions with unprecedented accuracy. There is therefore an immense 
interest in quantum resonances both in experimental and computational chemistry 
|Zar06l[5Y04llSSM+00| . 

We will show that the quantum normal form provides us with a very efficient 
algorithm for computing quantum resonances and also the corresponding resonance 
states. In our discussion of the classical reaction dynamics we could identify the ac- 
tivated complex with the centre manifold of the saddle-centre-- • • -centre equilibrium 
point, i.e. with an invariant subsystem with one degree of freedom less than the 
full system located between reactants and products (see Sec. 14. ip . As will discuss in 
detail in Sec. 16.31 the Heisenberg uncertainty relation excludes the existence of an 
invariant quantum subsystem. In fact, the quantum resonances will describe how a 
wavepacket initialised near the classically invariant subsystem will decay in time. 

Quantum resonance can be introduced in several ways. A common definition 
is based on the S-matrix. If one can extend the S-matrix analytically to complex 
energies, then the resonances are defined as its poles in the complex energy plane. 
We therefore could use the results of the previous section to determine the resonances 
from the quantum normal form. However, we will choose a different approach to 
introduce resonances which will make their dynamical meaning much more clear. 



6.1 Definition of quantum resonances 

We will define resonances as the poles of the resolvent operator. This is in line with 
the the convention in the mathematical literature (see, e.g., jZwo99j ). Let us recall 
the necessary notions. 

For an operator H : L^(]R'^) L^(R''), the resolvent set r{H) of H is defined as 
the set of ii^ G C such that H — E \s invertible. The spectrum of H is the complement 
of the resolvent set. For E £ r{H), the resolvent of H is defined as 

R{E) = {H- EY^ : l2(E^) ^ l2(R'^) . (6.1) 

If H is selfadjoint, then the spectrum of H is contained in R. The resolvent is thus 
defined at least for all E G C\1R,. The resolvent is related to the time evolution 
operator U{t) = exp(— -^t//) by Laplace transformation. For lm.E > 0, 



R{E) = ^ I et^*C/(t) dt , (6.2) 



and by Mellin transform 



Uit) = ^[ R{E)eT^^dE, (6.3) 



27ri 



where c > 0. The path of integration in the Mellin integral should be thought of 
as encircling the spectrum of H. Hence, if H has only isolated eigenvalues En then 
Cauchy's theorem gives 

C/(t)=^e-t*^"P„ (6.4) 



with the projectors 

Pn-=^ I R{E) dE, (6.5) 



27ri 



where the Cn are a closed paths encirchng only This is the usual spectral theorem 
which shows how eigenvalues and eigenfunctions (contained in the projectors Pn) 
determine the time evolution of a system with discrete spectrum. 

In the case where the spectrum of H is not discrete the sum over eigenvalues is 
replaced by an integral, and it becomes harder to read off properties of the time- 
evolution directly. Physically, a continuous spectrum corresponds to an open system 
like a scattering system where wavepackets can decay by spreading out to infinity. 
This will be described by resonances. 

Let us assume H has continuous spectrum. The resolvent R{E) is an analytic 
function of E for Imii^ > 0, and the resonances are defined as the poles of the 
meromorphic continuation of R{E) to the region ImE < 0. Since the operator 
H is selfadjoint on L^(R'^) and has continuous spectrum, there is no meromorphic 
continuation of R{E) as an operator from L^(IR'^) L'^(M.'^). Instead one looks for 
a continuation of R{E) as an operator 

R{E) : Ll^piR'') ^ LlM") , (6-6) 

where L'^g^p{W^) and L'^^^{li'^) denote the spaces of functions that are in L^(]R'^) and 
have compact support, or that locally are in L^{M!^), respectively. More directly, 
let ip,ij: £ L^Qj^p(IR'^), then quantum resonances are the poles of the meromorphic 
continuation of the matrix elements 

{f,R{E)i^) (6.7) 

from the region Im E > to Im E < 0. Assuming we have found such a meromorphic 
continuation with poles at En G C, n G M, Im En < 0, then we can use (16. 3p to get 

U{t)ij) = ^ [ iv, R{E)ij)eT.'^ dE . (6.8) 



27ri 



Im E=c 



Shifting the contour of integration and picking up the contribution from the poles 
gives us an expansion in terms of the resonances En 



with the projectors 



Pn--=^ ! R{E) dE, (6.10) 



27ri 



where C„ is a closed path encircling only the resonance En- This looks formally like 
()6.4p . but there are two important differences. Firstly, Im£'„ < which means that 
|g--^t£;„| _ gtimEn^ hence the terms in the sum are exponentially decreasing for 
t — > oo (since ImE'Tv < 0). Secondly, the projectors P„ are no longer orthogonal 
projectors in L'^(M!^). Futhermore, we can take the expansion only as far as the 
meromorphic continuation allows us to, and even if it extends to C, the resulting sum 
could be divergent. The range of the meromorphic continuation and the convergence 
properties of the sum can depend on if and ip (see |Zwo99j for a more detailed 
description). 

The relation ()6.9p reveals the dynamical meaning of the resonances. Resonance 
states are not stationary, and the reciprocal value of the imaginary part of the reso- 
nance energies determines their lifetime. 



6.2 Computation of resonances of the quantum normal 
form 



We now turn to explicit calculations and show how one can compute quantum reso- 
nances of the quantum normal form. 

6.2.1 Resonances of one-dimensional systems 

We start with the simplest one-dimensional example { d = I) and consider the oper- 
ator 

H = XI = X-(qdg + l), (6.11) 



where A > 0. For this operator the Schrodinger equation can be solved explicitly and 
the time evolution operator is given by 

U{t)iP{q) = e-^ip{e-^^q) . (6.12) 

This operator is of course unitary, i.e. it preserves the L^-norm. In time, the state 
U{t)'tp{q) spreads out at an exponential rate. If we look at the overlap of U{t)tp{q) 
with another localised state we expect an exponential decay, and this is exactly what 
the resonances describe. Let S C^(R), then 

{if, U{t)^) = y ip*{qme-^'q) dq (6.13) 

and if we insert for -0 its Taylor series 



n=0 

with |i?jv4.i(g)| < CN+i\(i\^^^ , then we obtain 

for t > Q. Inserting this equation into (j6.2p leads to the meromorphic continuation 
of R{E) to the domain ImE' > -h\{N + 1 + 1/2) with poles at 

En = -'ih\{n + 1/2) , n = 0,...,N . (6.16) 

These are the resonances of the operator H given in ()6.1ip . 
We can furthermore read off the projection operators 

Pni^iq) := -.^P^''H0)q\ (6.17) 
n! 

and a direct calculation shows that q^ is an eigenfunction with complex eigenvalue 
En = -mX{n + l/2), 

Hq^ = -m\{n + l/2)g" . (6.18) 

We now extend this analysis to the case of a Hamilton operator in quantum 
normal form for d = 1, i.e., H = K{I), where K is a polynomial or an analytic 
function in /. We will require furthermore the condition 

ImK(-ix) < , for x > . (6.19) 



By expanding K in a. power series we find 

fl-g" = -in(n + l/2))g" (6.20) 

and solving the Schrodinger equation yields U{t)q^ = exp[— i 
Hence, if il}{q) is analytic, we have 

oo ^ 

U{t)^{q) = —^lji^){{))e-W-ih{n+l/2))^n ^ 

and by condition (j6.19p we can use (j6.2p to see that the resonances are given by 

En = K{--\h{n + 1/2)) , n = 0,l,2, .... (6.22) 

This can be regarded as a kind of imaginary Bohr-Sommerfeld quantisation condition 
for the resonances. Moreover, we can formally write the resonance states ^n(<z) = 
as "complex" Lagrangian states, 

Mq) = = (sgn g)"|g|-V2+i/nM (6.23) 

with In = —ih{n + 1/2). This reveals the formal similarity of the resonance states to 
the scattering states (j5.5|) with the main difference being that in the case of resonances 
/ fulfils an imaginary Bohr-Sommerfeld quantisation condition while in the case of 
scattering the spectrum of / is continuous and real. With the states (j6.23|) we can 
associate the complex Lagrangian manifolds 

A0„ = {iq,p) = iqjn/p) : qeR} CRxiR. (6.24) 

6.2.2 Resonances of multi-dimensional quantum normal form 

Finally, we consider the case of a d-dimensional system in quantum normal form, i.e. 
let H = K{I, J2, . . . , Jd) and ^pn^ denote the n^^ harmonic oscillator eigenfunction 
(see (fOU]) l. For n = (m, . . . ,nrf) G INjj, set 

V'n(g) = ^1 Vn2(92) • ■■'fnAld) ■ (6-25) 

Then we have 

HtPn = K{- ih{ni + 1/2), ;i(n2 + 1/2), . . . , h{nd + 1/2)) Vn , (6.26) 

and if we assume Imi^(— ixi, X2 • • • , Xd) < for xi > and X2, • • • , in a neigh- 
bourhood of 0, we can conclude as before that the resonances of H are given by 

En = K{-ih{ni + 1/2), h{n2 + 1/2),..., h{nd + 1/2)) , n G . (6.27) 

To summarize, we have shown 

Theorem 4. Suppose H = K{i, J2, - ■ ■ , Jd) o.i'^d that K satisfies the condition 

ImE:(-ixi,X2,...,Xrf) < (6.28) 

for xi > and X2, ■ ■ ■ Xd in some neighbourhood of 0. Then the resonances in a 
neighbourhood of are given by 

En = K{-ih{ni + l/2),h{ni + l/2),...,h{nd + 1/2)) , n € < . (6.29) 

and the corresponding resonance eigenstates are 

i^niq) = ^1 Vnafe) • ■■fnaiqd) ■ (6.30) 

Following (j6.24p the resonance eigenstate can be interpreted as Lagrangian states 
associated with the complex Lagrangian manifolds 

A^„={{q,p) eR^'^ : Pi=InJqi, {pi + ql) = 2nkh , k = 2,...,d} . (6.31) 



6.3 Lifetime of the activated complex 

The geometric object in classical phase space associated with the activated complex 
is the centre manifold, a (2ci— 2)-dimensional invariant submanifold. As mentioned in 
Sec. l4.4l this submanifold can be considered as the phase space of a (d—l) DoFinvari- 
ant subsystem related to the supermolecule poised between reactants and products 
in the chemistry literature [Pec 76 1 IMar92j . This invariant subsystem is unstable, 
i.e. a trajectory with initial condition near but not in the subsystem will leave the 
neighbourhood of this subsystem. 

For the corresponding quantum system the Heisenberg uncertainty relation ex- 
cludes the existence of a quantum analogue of the classical invariant subsystem. This 
is because in normal form coordinates the invariant manifold is defined hy qi = pi =0 
and in quantum mechanics we have the uncertainty relation ApiAqi > h/2, i.e. pi 
and qi cannot be simultaneously. The closest one can get to a state which initially 
has (7i = pi = is a minimal uncertainty state which is a Gaussian of the form 

Mqi) = (^^^-^^ • (6-32) 



1piqi,---,qd) = ,_^^UA ^ ^ fn2iQ2) ■ ■ ■ fnaild) (6.33) 



In order to obtain a state which at time t = is localised on the centre manifold we 
choose 

1 

for some fixed quantum numbers n2, ■ ■ ■ ,nd £ INq, where (pn^. again denote the har- 
monic oscillator eigenfunctions. 

A suitable quantity for measuring the lifetime of such a state is the decay of the 
autocorrelation function 

\{i;,U{tm\ (6.34) 

We will compute the autocorrelation function for the case that the Hamiltonian is 
in quantum normal form. Inserting the expression (j6.33p for ip and expanding the 
Gaussian into a Taylor series gives 

^ °° I (—1)'^ 1 ^ 



_ 1 (-1)*' 1 f ISI 2k . -j;tH{-ih{2k+l/2),h{n2+l/2),-Mnd+im) 

(6.35) 

where we have used as well that q"^^ is a resonance state (j6.26p . The integral over qi 
gives f e'T^^qf^ dqi = T{k + l/2)(2/i)'^+^/2, and we thus find 

{lp,U{t)lp) = / 2 T{k + 1/2) ^ ^^fc^-^f J/(,m(2fc+l/2),a(n2+l/2),- ,h(na+l/2)) _ 

fc=0 ^' 

(6.36) 

The leading term in this sum for t — > cxd is given by the smallest resonance with 
fe = 0. Hence, 



H(^-ih/2,h{n2+l/2),--- ,S.(nd+l/2)) 



and this determines the maximal lifetime of a quantum state of the activated complex, 
i.e. a state initially localised on the invariant subsystem given by the centre manifold. 



For small h the quantum normal form is dominated by its quadratic part and 
that gives 



liin^21mH{-ih/2,h{n2 + 1/2),- •• ,h{nd + 1/2)) = -X (6.38) 
and therefore for small h 

K^,[/(t)V)|'~2e-*V (6.39) 

The quantum lifetime of the activated complex is in leading order for ^ ^ thus 
given by the reciprocal value of the classical Lyapunov exponent associated with the 
saddle equilibrium point. 

6.4 On the relation between the resonances of the Quan- 
tum Normal Form and the full system 

We have seen that the resonances of an operator in quantum normal form can be 
computed explicitly. They are obtained from the Bohr-Sommerfeld type quantization 
condition (I6.27p . In Sec. [3] we have shown how to approximate a Hamilton operator 
near an equilibrium point of the principal symbol by an operator in quantum normal 
form. We now want to discuss under which conditions this quantum normal form 
can be used to compute the resonances of the full Hamilton operator. This question 
has been studied in |KKOO] and we will mainly cite their results. 

One would expect that resonances of the full system are close to the one of 
the quantum normal form around an equilibrium point if that equilibrium point 
dominates the reaction, i.e., if it is the only equilibrium point at that energy, and all 
other trajectories come from infinity or can escape to infinity. This idea is formalized 
by using the trapped set of the classical Hamiltonian, whose definition we now recall. 

Let H{q,p) be a Hamilton function and the Hamiltonian flow generated by 
it. The trapped set at energy E is defined by 

TS^{H) := {(g,p) G x : H{q,p) = E, \ lim «>*j:^(g,p)| < oo} . (6.40) 

It consists of the trajectories which stay in some bounded region for t ±00. 

Theorem 5 ( [KKOOj ). Assume H satisfies the general conditions of IHS86^ and has 

a equilibrium point at zq with energy Eo and TS^°{H) = {zq}. Let K^^j^p he the 7V*^ 
order quantum normal form of H with respect to zq. Then the resonances of Op[H] 
in a h^ neighbourhood of Eq, 1 > 5 > 0, are close to the resonances of Kg^Jp. 

The conditions from |HS86j referred to above are conditions on H which ensure 
that the resonances can be defined by a complex deformation of phase space, a 
generalisation of the complex dilation method |Sim79[ IRei821 IMoi98j which we will 
use in Sec. [7] to compute numerically exact quantum resonances. For a more recent 
and more accessible presentation see [LBM02j . 

More explicitly, the main consequence of Theorem [5] is that for every n G INq, 
there is a resonance En G C of Op[//] with 



En = i^SU - + 1/2), /i(ni + 1/2), . . . , h{na + 1/2)) + 0{{\n\h)''+^) . (6.41) 

The quantum normal form thus provides an asymptotic expansion of the resonances 
for small h. If we want to have all resonances in a neighbourhood of Eq of radius , 
then we must go in n up to a size determined by ?i|n| ~ /i^ in which case the error 



term becomes of order h . Since we are interested in the first few resonances only 
we can take 6 = 1. 

We note that the resonances (j6.4ip coincide with the poles of the S-matrix which 
we computed in ()5.24p and ()5.32p . As can be seen from ()5.24p the poles of the S- 
matrix are simply given by the poles of the gamma function at non-positive integers. 

In cases when the trapped set is larger, e.g., when there are several equilbrium 
points at the same energy, the situation is more complicated and the structure of the 
set of resonance is no longer necessarily determined by the contributions from the 
individual equilibrium points. Instead one has to use the methods sketched in Sec. 
15.61 to construct a global S-matrix which will bring the global geometry into play. 

6.5 Distribution of the resonance states in phase space 

We now want to study the distribution of the resonance states in phase space in terms 
of Husimi functions. Like in the case of the scattering states in Sec. 15.51 the Husimi 
function of the resonance states (j6.25p is given by the product of the Husimi functions 
of harmonic oscillator eigenfunctions (pn,. and the Husimi function of (pmiQi) = Qi^ ■ 
We already discussed the Husimi function of the ipm. in Sec. 15.51 The computation 
of the Husimi function of the (pm is rather straightforward, and we obtain 

where is the nf^ Hermite polynomial. Therefore we have 

2 

e-P'/'^ . (6.43) 

Figure [12] shows contour plots of the Husimi functions of the first five resonance 
states. Due to the exponential damping in the direction of pi the Husimi functions 
H^^ are concentrated along pi = 0. Along pi = they increase in leading order in 
qi as 

/7^„(<?i,O)~-^0y q^ + Oiqr'). (6.44) 

It follows from ()6.43p that H^^_^ has rii zeroes located near the origin on qi = 0. 

For n = (ni, . . . , n^) € Wq the Husimi function of a multi-dimensional scattering 
wavefunction defined in (j6.25p is simply given by the product of the functions 
defined in (lOHIl and (fOH]) . i.e. 

H^pAq,?) = H^„^{qi,pi)H^^^{q2,P2) ■ ■ ■ H^^^{qd,Pd) ■ (6.45) 

From the distribution of the functions (I5.43P and ()6.43p it thus follows that the reso- 
nance states ipn are concentrated on the real projections of the complex Lagrangian 
manifolds A^^ in ()6.3ip 

{((?,p)gE2^ : pi=0, ipl + ql) = 2nkn, k = 2,...,d}. (6.46) 

7 Examples 

In the following we illustrate the classical and quantum reaction dynamics for concrete 
examples with one, two and three degrees of freedom. As we will see the reaction 
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Figure 12: Contour plots of the Husimi functions H^^ in the (gi,pi)-plane for ni = 0, . . . , 5. Red corre- 
sponds to low values; blue corresponds to high values. The spacing between the values of the contourlines 
is decreasing exponentially. [h = Q.l .) 



dynamics in systems with one or two degrees of freedom still has certain features that 
do not persist in the multidimensional case (of three or more degrees of freedom). We 
will use the classical normal form to realise the phase space structures that control 
classical reaction dynamics for these systems and compute the classical flux. Likewise 
we will use the quantum normal form to compute cumulative reaction probabilities 
and quantum resonances. We note that we implemented the procedures to compute 
the classical and quantum normal forms in the programming language C-I--I-. In our 
object-oriented implementation the number of degrees of freedom and the order of 
the normal form can be chosen arbitrarily. 



7.1 Example with 1 DoF 

The most frequently used systems to model one-dimensional reaction problems, like 
the paradigm hydrogen exchange reaction H2-I-H — > H -|- H2, are the parabolic barrier 
and the Eckart potential (see, e.g., |SM91[ [SY04]). The reason for choosing these 
model systems is that the reflection coefficient and the quantum resonances can be 
computed analytically for these systems. We have already seen that the quantum 
normal form computation of the reflection coefficient and the resonances is exact for 
a parabolic barrier. We therefore focus here on the Eckart barrier which provides a 
much more realistic model of reactions than the parabolic barrier. 
The Hamilton function for an Eckart barrier |Eck30j is given by 

H =p'/{2m) + V^{x), (7.1) 

where Ve is defined as 

/ N ^ exp((x xo)/a) exp((x -H xo)/a) 

V-E,[x) = A- -——+ B— — , (7.2) 

^ ' 1 + exp((a; + a;o)/a) (1 + exp((a; + a;o)/a))2 ^ ' 

with 

a;o = aln ^_^ . (7.3) 
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Figure 13: (a) Graph of the Eckart potential Ve defined in ()7.2p with parameters a = 1, B = 5 and 
different values of A. (b) Phase portaits for the Eckart potential with parameters a = 1, A = 0.5, B — 5 
and m = 1. The green and red lines mark the stable and unstable manifolds of the equilibrium point 



For B > A > the Eckart potential possesses a maximum which we shifted to x 
for convenience. The value of the potential at its maximum is 
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(7.4) 

oo (see 



The potential monotonically decreases to as x - 
Fig. [THb). For A = 0, the potential is symmetric. 

The Weyl quantisation of the Hamilton function ()7.ip gives the Hamilton operator 



Op[H] 



2m dx^ 



+ Ve. 



(7.5) 



The Hamilton function H in (17. ip is then the principal symbol of the Hamilton 
operator Op[ff]. 



7.1.1 Computation of the classical and quantum normal forms 

In order to compute the quantum normal form we can follow the calculation for one- 
dimensional potential barriers described in Sec. 13.51 Using the notation of Sec. 
the coefficients of the Taylor expansion to fourth order are 
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(7.6) 
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We refrain from giving the analytical expressions for the higher order terms as the 
actual computation of the classical and quantum normal form implemented in our 
C++ program is carried out numerically. However, we used the coefficients above 
together with Equation ()3.143p to check the numerically computed 4**^ order quantum 
normal form. To give the reader the opportunity to verify our results we list in Tab. [Ill 
in Sec. [B] of the appendix the coefficients of the symbol of the 10*^ order quantum 
normal form of the Eckart barrier with parameters a = 1, B = 5, A = 1/2 and 
m = 1. The classical normal form can be obtained from the symbol by discarding all 
terms that involve a factor h. 



7.1.2 Classical reaction dynamics 



Since the energy surface of a 1 DoFsystem is one-dimensional, the classical reaction 
dynamics of 1 DoFsystems is trivial. The question of whether a trajectory is reactive 
or nonreactive is determined by the energy alone, i.e., in the case of the Eckart barrier 
trajectories are forward or backward reactive if they have energy E > Ve{0), and they 
are nonreactive localised in reactants or products if ii^ < 1^(0) (see Fig. 113b). Fixing 
an energy E > Ve{0) one can choose any point G R to define a dividing 'surface' 
on the energy surface according to {{x,px) ■ x = xjs, H{x,px) = E} = {{x,px) = 
(xds, ±\/2m{E — Ve{x))}. This dividing 'surface' consists of two points which have 
Px > and Px < and are crossed by all forward reactive trajectories and backward 
reactive trajectories, respectively. In fact the two points can be considered to form 
a zero-dimensional sphere, with each point forming by itself a zero-dimensional 
ball, 5°. Note that many of the other phase space structures that we discussed in 
Sec. m do not make sense for the case of d = 1 degrees of freedom. Moreover, the case 
of one degrees of freedom is special because it is the only case for which the location 
of the dividing surface is not important. 

Note that the formalism to compute the classical flux f{E) developed in Sec. 14.41 
does not apply either to the case d = 1. Still it useful to view the classical flux to 
be given by the step function f{E) = @{E — Ve{0)), i.e., classically, we have full 
transmission for E > Ve(0) and full reflection for E < Ve(0). 



7.1.3 Quantum reaction dynamics 



The effect of quantum mechanical tunneling makes the quantum reaction dynamics 
even of 1 DoFsystems more complicated than the corresponding classical reaction 
dynamics. The quantum mechanically exact transmission coefficient Texact can be 
computed analytically for the Eckart potential [Eck30| . One finds 



where 
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Note that Texact (-E') — > when the energy E approaches the limiting value A of the 
potential from above. Figure shows the graph of Texact (-E') versus the energy E. 
Following Sec. 15.31 we can compute the transmission coefficient from the A^**^ order 



quantum normal form -K^qnf according to 
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where I^^\E) is obtained from inverting the equation 
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(7.10) 



(7.11) 



We illustrate the high quality of the quantum normal form computation of the 
transmission coefficient in Fig. [T5k which shows the difference between Texact and 
Tq^p for different orders, N , of the quantum normal form. Though the quantum 
normal form expansion is not expected to converge, the difference between Texact and 
Tq^p decreases as increases to the maximum value of 10 at which we stopped the 
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Figure 14: Exact transmission coefficient Tcxa.ct{E) (top panel) and resonances in the complex energy 
plane (bottom panel) for the Eckart potential. The parameters a, B, C and m are the same as in Fig. 113b . 
and h = 0.1. 



quantum normal form computation. In fact, the difference decreases from the order 
of 1 percent for the 2^^*^ order quantum normal form to the order of 10~^^ for the 
■j^Qth Qp^jgj. quantum normal form. 

We can also compute the quantum mechanically exact resonances analytically. 
They are given by the poles of the transmission coefficient (j7.8|) . We find 

-C/exact,n — (-^ TZ r, TTT5 ) n — U,i,^,.... I'-J-^j 

(0 - i(n + 

We illustrate the location of the quantum resonances in the complex energy plane 
in the bottom panel of Fig. [TH Following Sec. 16.2.11 we can compute the resonances 
from the A^**^ order quantum normal form according to 

4NF,n = ^QNF(-i^(^ + 1/2)) , n = 0, 1, 2, . . . . (7.13) 
For the 2"^^^ order quantum normal form this reduces to 

^QNF,n = ^E(0)-iA;i(n + i), n = 0,1,2,... . (7.14) 

As mentioned earlier the 2°^^ order quantum normal form resonances would be exact 
for a parabolic potential barrier. For comparison we also show the location of these 
resonances in the complex energy plane in the bottom panel of Fig. [TH Note that 
the 2°"^ order resonances have a constant real part. The "bending" of the series 
of exact resonances in Fig. [T3] is a consequence of the nonlinearity of the Eckart 
potential. The quantum normal form is able to describe this effect very accurately. 
The approximation of the exact resonances by the 4*^ order quantum normal is 
already so good that the error is no longer visible on the scale of Fig. [TH We therefore 
show the differences between the exact and quantum normal form resonances for 
different orders of the quantum normal form in a separate graph in Fig. [TSb . Again, 
up to the maximal order shown, the accuracy of the quantum normal form increases 
with the order. As to be expected, for a fixed order of the quantum normal form A^, 
the error of the quantum normal form increases with the quantum number n. Note 
that the sequence of resonances is localised in the complex energy plane in Fig. [T3] in 
such a way that the real part of the resonance closest to the real axis coincides with 
the position of the (smooth) step of the transmission coefficient on the (real) energy 
axis. 



a) b) 




E 



Figure 15: (a) Error for the transmission coefficient of the Eckart potential computed from quantum 
normal forms of different orders N. (b) Errors for the the resonances of the Eckart potential computed 
from quantum normal forms of different orders as a function of the quantum number n. The parameters 
for the Eckart potential are the same as in Fig. [T31 



7.2 Example with 2 DoF 

We now illustrate the quantum normal form computation for a 2 DoFmodel system 
which consists of an Eckart barrier in the x-direction that is coupled to a Morse 
oscillator in the y-direction. A Morse oscillator is a typical model for a chemical 
bond. The Hamilton function is 

H = i^{pI+pI) + VEix) + Vuiy) + eH, , (7.15) 

where Ve is the Eckart potential from ()7.2p and Vm is the Morse potential 

Vuiy) = De(exp(-2aMy) - 2exp(-aMy)) (7.16) 

with positive valued parameters Dg (the dissociation energy) and om (see Fig. 116b,). 
For the coupling term He we choose a so called kinetic coupling (see, e.g., jHel95| ) 

Hc=PxPy (7.17) 

The strength of the coupling is controlled by the parameter e in (I7.15p . The vec- 
tor field corresponding to the Hamilton function (j7.15p has an equilibrium point at 
{x,y,Px,Py) = 0. For |e| sufficiently small (for given parameters of the Eckart and 
Morse potentials), the equilibrium point is of saddle-centre stability type. Contours 
of the Eckart-Morse potential V{x, y) = Ve{x) + VM{y) are shown in Fig. [TCb . These 
indicate the bottleneck-type structure of the energy surfaces with energies slightly 
above the energy of the saddle-centre equilibrium point. Note that the relation be- 
tween the saddle of the potential V{x,y) = V-e,{x) + Vyiiy) at (x,y) = and the 
equilibrium point of Hamilton's equations at {x,y,Px,Py) = is complicated by the 
kinetic coupling in (I7.15p . 

The Weyl quantisation of the Hamilton function H in (|7.15p gives the operator 



The Hamilton function ()7.15p is the principal symbol of the operator Op[-ff]. 



Figure 16: (a) Morse potential Vuiy) — (exp(— 2aM2/) — 2 exp(— 2aM y))- The potential approaches 
for 2/ ^ cxD. The parameter — Vivi(oo) — Vm(0) gives the depth of the well potential well while 
om determines the width of the well, (b) Contours of the Eckart-Morse potential Ve{x) + Vuiy)- Red 
correspond to small values of the potential; blue corresponds to large values. The parameters for the Eckart 
potential are the same as in Fig.[T4l The parameters for the Morse potential are De = I and au = 1- 



7.2.1 Computation of the classical and quantum normal forms 



Since the equilibrium point is already at the origin of the coordinate system we 
can skip the first step in the classical and quantum normal form transformation se- 
quences ()2.27p and (j3.65p , and start with the second step which consists of simplifying 
the quadratic part of the Hamilton function or symbol, respectively. To this end we 
follow Sec. 12.31 and compute the matrix JD^H associated with the linearisation of 
Hamilton's equations about {x,y,px,Py) = 0. This gives 
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(7.19) 



where Ae is defined as in (17.60 and 
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is the linear frequency of the Morse oscillator. The matrix in (17.190 has eigenvalues 
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where as mentioned above, for given parameters of the Eckart and Morse potentials 
and |e| sufficiently small, the eigenvalues ei and 63 (and hence A) are real, and 62 
and 64 are purely imaginary (and hence lo is real). For e ^ 0, A and uj converge to 
Ae and ujm, respectively. The corresponding eigenvectors are 



Vk 



(/O 0\ O 0/0 0\ 000\ 7^ 

ek[ei^+uj ),ernX-^ek,mX-^[ek + uj ),-ern X-p^uj ) , A; = 1,2,3,4. (7.25) 



Following Sec. 12.31 we obtain a real linear symplectic change of coordinates by using 
the Vk to define the columns of a matrix M according to 



M = (cit;i,C2Ret;2,ciW3,C2lmt;2) 



(7.26) 



with the coefficients ci and C2 defined as 
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(7.27) 



Now set 



(7.28) 



Then the Hamilton function (I7.15P becomes 



H = V{Q) + \qipi + -{ql+pf, + . . . 



(7.29) 



where the neglected terms are of order greater than 2. The constant term is 



The truncation of (j7.29p at order 2 is the symbol of the 2^^*^ order quantum normal 
form of (TTTHD . 

The classical and quantum normal forms are then computed from the algorithm 
described in Sections 12.31 and respectively. For the parameters a = 1, i? = 5, 
A = 1/2 for the Eckart potential and Z^e = 1 and om = 1 for the Morse potential, 
e = 0.3 for the coupling strength, and m = 1, we list the coefficients of the symbol 
of the 10**^ order quantum normal form in Tab. IIIII of the Appendix. The classical 
normal form can be obtained from the symbol by neglecting all terms that involve a 
factor h. 

7.2.2 Classical reaction dynamics 

For a 2 DoFsystem the NHIM is a one-dimensional sphere, S^, i.e., a periodic orbit. 
This is the Lyapunov periodic orbit associated with the saddle point. As discussed in 
the introduction, for 2 DoFsystems with time-reversal symmetry, the periodic orbit 
can be used to define a dividing surface without recrossing - the so called periodic 
orbit dividing surface - from the projection of the periodic orbit to configuration space 
|PM73[|PP78] . Note that, as mentioned earlier, such a construction in configuration 
space does not work for systems with 3 or more DoF [ WW04j . 

The NHIM has stable and unstable manifolds with the structure of cylinders or 
'tubes', X R. They inclose the forward and backward reactive trajectories as 
discussed in detail in, e.g., |WBWn5bl TWEWOScj . The flux is given by the action of 
the periodic orbit |WW04| . In the uncoupled case (e = 0) the periodic orbit (p.o.) is 
contained in the (7/,py)-plane and its action can be computed analytically. One finds 



^(0) = yE(0)+yM(0) 
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(7.31) 



for -De + Ve{0) < E < Ve{0) and f{E) = (no classical transmission) for E < 

-De + VE{0). 




Figure 17: The top panel shows the cumulative reaction probabilities iVcxact(-E') (oscillatory curve) and 
7Vweyi(^') for the Eckart-Morse potential defined in the text with e = 0. The bottom panel shows the 
(numerically) exact resonances computed from the complex dilation method in the complex energy plane. 
Circles mark resonances for the uncoupled case e = and crosses mark resonances for the strongly coupled 
case e = 0.3. The parameters for the potential are the same as in Fig. \Wi Again we choose m = 1 and 
h = 0.1. 



7.2.3 Quantum reaction dynamics 

For the uncoupled case we can compute the cumulative transmission probability 
analytically. We have 

^exact{E) = ^Eckart; exact (-^^ — -EMorse;n2) ) (7.32) 

where TEckart; exact denotes the transmission coefficient for the Eckart barrier given 
in (|7.8|) and -E'Morse;n2 energy levels of a one-dimensional Morse oscillator 



£^Morse;n2 = U2 + - - , 712 = 0,1,2,.... (7.33) 



2m V 2 auh 

The graph of A'^exact in the top panel of Fig. [T7] shows that A'^exact is "quantised" , 
i.e., it increases in integer steps each time a new transition channel opens. The 
opening of a Morse oscillators mode (7x2) as a transition channel can be defined 
as the energy where Tsckart; exact - -E'Morse;n2) = 1/2- The quantisation of the 
cumulative reaction probability has been observed experimentally, e.g., in molecular 
isomerisation experiments |LM93j and also in ballistic electron transport problems 
in semiconductor nanostructures where the analogous effect leads to a quantised 
conductance |vWvHB+88| IWTN+88] . As mentioned in Sections [Ol and [5^ the 



quantity 

NweyliE) = f{E)/{2TTh) (7.34) 

can be interpreted as the mean number of open transmission channels at energy E. 
This is illustrated in the top panel of Fig. [T7] which shows A'^vVeyi together with N exact- 
Note the nonlinear increase of A'weyi(-E') with E which is an indication of the strong 
anharmonicity of the Morse oscillator. 



In order to compute the cumulative reaction probability from the quantum normal 
form we follow the procedure described in Sec. 15.41 We get 

-1 

(7.35) 

where Inl\E) is obtained from inverting 

^QNF(^nf (i^),?i('^2 + l/2)) =i?, 712 = 0,1,2,.... (7.36) 

The high quality of the quantum normal form computation of the cumulative reaction 
probability is illustrated in Fig. [T8b which shows | A^qnf {E) — iVexact {E) \ versus the 
energy E for quantum normal forms with = 2 to = 10. Like in the 1 DoFexample 
in Sec. 17. II we find that up to the orders shown, the accuracy of the quantum normal 
form increases with the order of the quantum normal form. The error is of order 
-j^Q-io £qj. j^Qth Qpjjgj. quantum normal form. 

For the coupled case e 7^ we also make a comparison of the quantum mechani- 
cally exact resonances and the resonances computed from the quantum normal form. 
The exact resonances cannot be computed analytically for the coupled case. To get 
them numerically we use the complex dilation method [Sim79l IRei82l IMoi98j whose 
implementation for the present system we describe in Sec. Oof the Appendix. The 
bottom panel in Fig. [T7] shows the (numerically) exact resonances for the uncoupled 
case and the strongly coupled case e = 0.3. In both cases the resonances form a dis- 
torted lattice in the complex energy plane. The quantum normal form computation 
of the resonances is given by 

^SWni,n2)=^QNFH^K + V2),M^2 + l/2)), ^ , n2 = 0, 1 , 2, . . . . (7.37) 

One of the benefits of the quantum normal form is that it leads to an assignment 
of the resonance lattice by quantum numbers. The quantum number ni labels the 
resonances in vertical direction, and the quantum number n2 labels the resonances in 
horizontal direction. Each vertical string of resonances (i.e., sequence of resonances 
for fixed 712) gives rise to one quantisation step of the cumulative reaction probability. 
Note that an assigment of the resonances is very difficult to obtain only from the exact 
quantum computation. Figure [T8b illustrates the high accuracy of the quantum 
normal form computation for a selection of resonances. 
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7.3 Example with 3 DoF 

Our final example is a 3 DoFmodel system consisting of an Eckart barrier in the x- 
direction that is coupled to Morse oscillators in the y-direction and in the z-direction. 
The Hamilton function is 

H=^{pI+pI+ pI) + Ve{x) + VM;2{y) + Vm;3{z) + eH, , (7.38) 

where Ve is the Eckart potential from ()7.2|) and VM;k, k = 2,3, are Morse potentials 
of the form (j7.16p with parameters D^-k and aj^f./^, k = 2,3, respectively. For He we 
choose the mutual kinetic coupling 

He =PxPy + Px Pz +PyPz- (7.39) 



The strength of the coupling is again controlled by the parameter e in (|7.38|) . The vec- 
tor field generated by the Hamilton function has an equilibrium point at {x,y, z,px,Py,Pz) = 
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Figure 18: (a) Errors for the cumulative reaction probability in the top panel of Fig. [T7] for different 
orders N of the quantum normal form, (b) Difference |-Eq^f — -Ecxactl for a selection of resonances with 
quantum numbers (ni, ^2) for the resonances shown in the bottom panel of Fig. [17] for the coupled case 
e = 0.3. 



0. For |e| sufficiently small (for given parameters of the Eckart and Morse potentials), 
the equilibrium point is of saddle-centre-centre stability type. Figure [19] shows con- 
tours of the potential V{x,y,z) = Ve{x) + Vm;2(2/) + Vm;3(^) which, for energies 
slightly above the saddle-centre-centre equilibrium point, indicate the bottleneck- 
type structure of the corresponding energy surfaces. 

The Weyl quantisation of the Hamilton function H in (|7.38p gives the operator 

^2 / q2 q2 q2 \ ^ / g2 q2 q2 \ 

^ 2m \dx^ ~'~ ~'~ dz^)''^ '^'^ \dxdy ~'~ dxdz ~'~ dydz J 

(7.40) 

The Hamilton function (|7.15p is the principal symbol of the operator Op[i7]. 



7.3.1 Computation of the classical and quantum normal forms 



Like in Sec. I7.2] the equilibrium point is again already at the origin of the coordinate 
system. For the computation of the classical and quantum normal forms we therefore 
again start with the second step in the sequences (|2.27|) and (|3.65p . respectively. 
Following again Sec. 12.31 we compute the Hamiltonian matrix associated with the 
linearisation of Hamilton's equations about the equilibrium point (x, y, z,px,Py,Pz) = 
0. This gives 
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where Ae is defined in (I7.6|) and 
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k — 2,3. 



(7.41) 



m 



(7.42) 



are the linear frequencies of the Morse oscillators. The matrix JD^H{0) has six 
eigenvalues, one pair of real eigenvalues of opposite signs and two pairs of imaginary 




Figure 19: Contours Ve{x) + VMaiv) + Vj\/;3(z) =const. of the Eckart-Morse-Morse potential. The 
parameters for the Eckart potential are the same as in Fig. [TH The parameters for the Morse potentials 
are D^a = ana = aM-.a = 1 and Dg^a = 2/3. 



eigenvalues with opposite signs. We label them according to 

ei = A , 64 = -A , 62 = ia;2 , 65 = -\uj2 , 63 = iws , 66 = -1^3 , 



(7.43) 



where A, U2 and uj^ are real positive constants that converge to Ae and the linear 
frequencies ljm;2 and wm;3, respectively, when e ^ . We assume that the parameters 
De-k and aM;fc, k = 2,3, are chosen such that uj2 and 0^3 are linearly independent over 
li. Let us again denote the corresponding eigenvectors hy v^, k = 1, . . . ,6. In order 
to define a real linear symplectic change of coordinates we use the eigenvectors Vk to 
define the columns of a matrix M according to 



M = {ciVi,C2 Re 172, 63 Re?;3,cif4,C2 Imt>2, C3 Imf3) 
with the coefficients ci, C2 and C3 defined as 



c-^ ■.= {vi,Jvi), C2 := (Ref2, JImt;2) 



(Re ^3, Jim ^3) 



(7.44) 



(7.45) 



We choose the eigenvectors vi and ^3 such that {vi, Jv^) is positive (if {vi, JV2) < 
then multiply V2 by -1). As mentioned in Sec. [ 
(|7.45|) are automatically positive and the matrix M is symplectic. For 



the coefficients ^ and Cg ^ in 



(gi,92,g3,Pi,P2,P3)^ = M ^{x,y,z,p^,py,pzf 
the Hamilton function (|7.38p becomes 



where the neglected terms are of order greater than 2. The constant term is 
V{0) = Ve{0) + Vm-2{0) + Vm;3(0) = ^^^^ - De-2 - D,.s . 



(7.46) 
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Figure 20: The NHIM projected into configuration space. The energy is 0.1 above the energy of the 
saddle-centre-centre equilibrium point. 

The truncation of (|7.47|) at order 2 is the symbol of the 2°*^ order quantum nor- 
mal form of (I7.38p . The higher order classical and quantum normal forms are then 
computed from the algorithm described in Sections 12.31 and 13. 3i For the parameters 
a = 1, B = 5, A = 1/2 for the Eckart potential and Dg-i = 1, De-2 = 3/2 and 
o-M;i = 0'M]2 = 1 for the Morse potential, e = 0.3 for the coupling strength, and 
m = 1, we list the coefficients of the symbol of the lO*'^ order quantum normal form 
in Tab. IIVI of the appendix. The classical normal form is obtained from discarding 
terms involving a factor h. 

7.3.2 Classical reaction dynamics 

The NHIM is a three-dimensional sphere, S^. In Fig. [20] we show the NHIM with the 
energy 0.1 above the energy of the saddle-centre-centre equilibrium point projected 
into configuration space, with the equipotential at the same energy for reference. 
Note that the projection of the NHIM to configuration space is a three-dimensional 
object. This can be viewed as an indication that the construction of a (in this case 
two-dimensional) dividing surface 'in configuration space' without recrossing is not 
possible for a system with 3 (or more) DoFsince, as explained in detail in [WW04] . 
a dividing surface without recrossing needs to contain the NHIM (as its equator). 

The NHIM's stable and unstable manifolds have the structure of spherical cylin- 
ders, X R. In Fig. [21] we show projections into configuration space of local pieces 
of the backward branch of the stable manifold of the NHIM, the forward branch of 
the stable manifold of the NHIM, the backward branch of the unstable manifold of 
the NHIM, and the forward branch of the unstable manifold of the NHIM. Due to 
the time-reversal symmetry of the system the stable and unstable manifolds project 
onto each other in configuration space. The stable and unstable manifolds enclose 
the forward and backward reactive trajectories as discussed in Sec. [U 

The NHIM is foliated by invariant 2-tori. According to Sec. 14.41 the classical flux 
for an energy E is given by 

f{E) = {2TifV{E) , (7.49) 



Figure 21: The stable and unstable manifolds of the NHIM projected into configuration space. Due to 
the time-reversal symmetry, these manifolds project onto each other in configuration space. The two colors 
represent the forward and backward branches of the manifolds, and they are "joined" at the NHIM. The 
energy is 0.1 above the energy of the saddle-centre-centre equilibrium point. 
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Figure 22: Energy contours in the plane of the Morse oscillator actions (J2, J3). The morse oscillators 
have energies Ey and Ez such that Ve{Q) + Ey + Ez = E with Ve(0) being the height of the one-dimensional 
Eckart barrier (see text). The parameters for the potential are the same as in Fig. [191 The mass m is 1. 




V{E) = I j^^E-VE{0)-Ey)^^dEy (7.51) 



where V{E) is the area enclosed by the energy contour in the plane of the corre- 
sponding action variables J2 and J3. In the uncoupled case the 2-tori are given 
by the cartesian products of two circles that are contained in the (y,py)-plane and 
{z,pz)-plane, respectively. The corresponding action variables J2 and J3 can be easily 
computed in this case. Let Ey and Ez be the energies contained in these two DoF. 
Then 

J2iEy) = ^(f pydy = — (^2mDe;2 - V-'^^Ey) > -De<Ey<0, (7.50) 

Jp.o. ■ «2 

and similarly for J^{Ez). The NHIM has energy E = V-e,{Q) + Ey + Ez, where 
Ve(0) = [A + B)'^/{4B) is the height of the one-dimensional Eckart barrier. Fig. 
shows some energy contours in the (J2, J3)-plane. The fact that the energy contours 
are no straight lines is an indication of the strong nonlinearity of the Morse oscillators 
for the energies shown. For an energy Ve{0) — De;2 — Dg-^ < E < Ve{0) — De-^, the 
inclosed area is given by 

E+D,.3-Ve(0) dJo(E ) 

= ^^^^^(a/^ - JVe{0) - D^,, - E) 
0.2 ^ V 

nrn 

{g{E - Ve{Q) + De;3) - g{-De;2)) , (7.52) 

02 as 

where 

(7.53) 

For E < Ve{0) — De-2 — De-^ the classical flux is zero. The graph of N^^Qyl{E) = 
f {E) / {2T:lif is shown in the top panel of Fig. [23l 

7.3.3 Quantum reaction dynamics 

In the uncoupled case the exact cumulative reaction probability A'^exact can be com- 
puted analytically. We have 

^exact(-E) = ^ ^Eckart; exact 

{E — EMorse;2,n2 ) , (7.54) 

"2, "3 

where TEckart; exact denotes the transmission coefficient for the Eckart barrier given 
in ()7.8p and £'Morse;k,nfc , k = 2,3, are the energy levels of the one-dimensional Morse 
oscillators, 



£^Morse;k,nfe = ^ h^fc + TT r~ ' rifc = 0, 1, 2, . . . . (7.55) 

2m \ 2 aM;kri 



The graph of A^exact gives the oscillatory curve shown in the top panel of Fig. | 

For the quantum normal form computation of the cumulative reaction probability 
we get 

(7.56) 



112, ns 



where I^^J ^^^{E) is obtained from inverting 

KQMlinlns)i^)^H^2 + 1/2), h{ns + 1/2)) = E , 712,713 = 0,1,2,... . (7.57) 

The high quahty of the quantum normal form approximation of the cumulative reac- 
tion probability is illustrated in Fig. [2lb which shows \N^p{E) — NcxactiE)\ versus 
the energy E. 

For the coupled case e 7^ we again make a comparison of the exact resonances 
and the resonances computed from the quantum normal form. We again compute 
the (numerically) exact resonances from the complex dilation method whose imple- 
mentation is described in Sec. Oof the Appendix. The bottom panel in Fig. 1231 shows 
the exact resonances for the uncoupled case and the strongly coupled case e = 0.3. 
In both cases the resonances now form a superposition of distorted lattices. The 
quantum normal form computation of the resonances 

^QNF (ni n, n,) = ^QNF(-i^(«l + V2), + 1/2), ^7(773 + 1/2)) , 77i , 772 , 773 G Mq, 

(7.58) 

allows one to organise the resonance by quantum numbers. The quantum numbers 
77i label the resonances in vertical direction, and the pairs of Morse oscillator mode 
quantum numbers (772,773) label the resonances in horizontal direction. Here each 
vertical string of resonances (i.e., sequence of resonances for fixed (772,713)) gives rise 
to one step of the cumulative reaction probability. In the top panel of Fig. [23] we 
mark the energies at which a mode (772,773) opens as a transmission channel. These 
energies are defined in the same way as in Sec. 17. 2i 

Since the density of the resonances in the complex energy plane is higher for the 
3 DoFcase than it is in the 2 DoFcase the quantisation of the cumulative reaction 
probability is more "washed out" . Again note that an assigment of the resonances is 
very difficult to obtain only from the exact quantum computation. The resonances 
computed from the quantum normal form are again of a very accuracy as shown for 
a selection of resonances with quantum numbers (771,772,773) in Fig. 124b. 



8 Conclusions and Outlook 

In this paper we have developed a phase space version of Wigner's dynamical tran- 
sition state theory for both classical and quantum systems. In the setting of Hamil- 
tonian classical mechanics, reaction type dynamics is induced by the presence of a 
saddle-centre-- •• -centre equilibrium point ('saddle' for short). For a fixed energy 
slightly above the energy of the saddle, the energy surface has a wide-narrow-wide 
structure in the neighbourhood of the saddle. Trajectories must pass through this 
bottleneck in order to evolve from reactants to products. We provided a detailed 
study of the phase space structures which for such an energy, exist near the sad- 
dle and control the dynamics in the neighbourhood of the saddle. In particular we 
showed the existence of a dividing surface which is free of local recrossings, i.e. it 
has the property that all trajectories extending from reactants to products (or vice 
versa) intersect this dividing surface exactly once without leaving a neighbourhood of 
the saddle and nonreactive trajectories which enter the neighbourhood from the side 
of reactants (resp. products) and exit the neighbourhood back to reactants (resp. 
products) do not intersect the dividing surface. This dividing surface minimizes the 
directional flux in the sense that a (generic) deformation of the dividing surface leads 
to an increase of the directional fiux through the dividing surface. Such a dividing 
surface is a prerequisite for the computation of reaction rates from the directional 
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Figure 23: The top panel shows the cumulative reaction probabilities iVoxact(-E') (oscillatory curve) and 
Nwcy\{E) (smooth curve) for the Eckart-Morse-Morse potential defined in the text with e = 0. It also 
shows the quantum numbers (^2,^3) of the Morse oscillators that contribute to the quantization steps. 
The bottom panel shows the resonances in the complex energy plane marked by circles for the uncoupled 
case e = and by crosses for the strongly coupled case e = 0.3. The parameters for the Eckart potential 
are the same as in Fig. [TH The parameters for the Morse potential are De-2 — 1, D^-^ = 3/2, aM:2 = 1 
and aM;3 = 1- Again we choose m = 1 and h = 0.1. 
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Figure 24: (a) Errors for the cumulative reaction probability in the top panel Fig. [^3] for different orders 
N of the quantum normal form, (b) Errors |i?Q^p — ^-cxactl for a selection of resonances with quantum 
numbers {ni, n2, 113) for the coupled case e = 0.3 in the bottom panel of Fig. 1231 



flux and its construction for multi-degree-of-freedom systems was considered a major 
problem in transition state theory. We showed that the existence of such a dividing 
surface is related to the presence of a normally hyperbolic invariant manifold (NHIM) 
which exists near the saddle. The NHIM has the structure of a sphere of two dimen- 
sions less than the energy surface. It can be considered to form the equator of the 
dividing surface which itself is a sphere of one dimension less than the energy surface. 
This way the NHIM divides the dividing surface into two hemispheres of which one 
is intersected by all trajectories evolving from reactants to products and the other is 
crossed by all trajectories evolving from products to reactants. 

The NHIM is the mathematical manifestation of what is referred to as activated 
complex in the chemistry literature. In fact, the NHIM, which is the intersection 
of the centre manifold of the saddle with the energy surface of the (full) system, 
can itself be viewed as the energy surface of an unstable invariant subsystem (the 
subsystem given by the centre manifold) . This subsystem has one degree of freedom 
less than the full system and as a kind of super molecule is poised between reactants 
and products. The theoretical background presented in this paper thus shows that 
the activated complex is not merely an heuristic concept utilised by transition state 
theory, but a geometric object of a precise significance for the dynamics. In particular, 
the NHIM has stable and unstable manifolds which have sufficient dimensionality to 
act as separatrices. They form the phase space conduits for reactions in the sense 
that they enclose the reactive volumes (which consist of trajectories evolving from 
reactants to products or vice versa) and separate them from the nonreactive volumes 
(which consist of nonreactive trajectories). They have the structure of spherical 
cylinders (i.e., cylinders where the base is a sphere). We discussed how the centre 
lines of the reactive volumes enclosed by these spherical cylinders naturally lead to 
the definition of a reaction path, kind of guiding trajectories about which 

other reactive trajectories rotate in phase space (observed as an oscillation when 
projected to configuration space) in a well defined manner. In contrast to the usual, 
often heuristic definitions of a reaction path, the reaction path presented in this paper 
incorporates the full dynamics in a mathematically precise way. 

We showed that all the phase space structures mentioned above can be realised 
through an efficient algorithm based on a standard Poincare-Birkhoff normal form. 
This algorithm allows one to transform the Hamilton function which describes the 
classical reaction dynamics to a simpler ('normal') form to any order of its Taylor 
expansion about the saddle point through a succession of symplectic transformations. 
In several examples we showed that the normal form computation truncated at a 
suitable order leads to a very accurate description of the dynamics near the saddle. 
In the generic situation where there are no resonances between the linear frequencies 
associated with the centre direction of the saddle the normal form is integrable and 
explains the regularity of the motion near the saddle which has been discovered in the 
chemistry literature |HB93[ IKB991 IMil77j . The integrability leads to a foliation of the 
neighbourhood of the saddle by invariant Lagrangian manifolds. These Lagrangian 
manifolds have the structure of toroidal cylinders, i.e. cylinders where the base is 
formed by a torus. 

We showed that similar to the unfolding of the classical dynamics in the neigh- 
bourhood of a saddle point we can obtain an unfolding of the corresponding quantum 
dynamics. We therefore reviewed some basic tools from the theory of micro local anal- 
ysis which allow one to study properties of quantum operators in a region of interest 
in the phase space of the corresponding classical system. The main idea is to use the 
Weyl calculus to relate Hamilton operators to phase space functions (symbols) and 
vice versa. This way one can extract properties of a Hamilton operator resulting from 



some classical phase space region by studying its symbol restricted to (or 'localised 
at') this phase space region. In the case of reaction dynamics the region of interest 
is the neighbourhood of a saddle point. We showed that in the neighbourhood of a 
saddle the Hamilton operator can be transformed to a simple form - the quantum 
normal form - to any order of the Taylor expansion of its symbol about the saddle 
by conjugating the Hamilton operator by a succession of suitable unitary transfor- 
mations. We showed that the quantum normal form computation can be cast into an 
explicit algorithm based on the Weyl calculus. This algorithm consists of two parts 
of which the first part takes place on the level of the symbols and is therefore very 
similar in nature to the classical normal form computation. The main difference is 
that the Poisson bracket involved in the symplectic transformation in the classical 
case is replaced by the Moyal bracket. In the second part of the quantum normal form 
algorithm the symbols are quantised to obtain the corresponding quantum operators. 
For this part we also developed an explicit algorithm. 

Through applications to several examples we illustrated the efficiency of the quan- 
tum normal form algorithm for computing quantum reaction quantities like the cu- 
mulative reaction probability and quantum resonances. The cumulative reaction rate 
is the quantum analogue of the classical flux. Quantum resonances describe the de- 
cay of wavepackets initialised on the centre manifold. In fact, quantum mechanically, 
the Heisenberg uncertainty principle excludes the existence of an invariant subsystem 
representing the activated complex analogously to the classical case. So quantum me- 
chanically a state initially localised on the centre manifold is unstable and will spread 
out. The quantum resonances describe the lifetimes of such states. We showed that 
these resonances are also related to the stepwise increase ("quantisation") of the 
cumulative reaction probability as a function of energy. The dependence of the cu- 
mulative reaction probability on the energy and also the resonances are viewed as the 
quantum signatures of the activated complex, and there is a huge experimental in- 
terest in these quantities |SY04j . In fact, recent advances in spectroscopic techniques 
allow one to study quantum scattering with unprecedented detail (see, e.g., |Zar06j ). 
We hope that the results presented and the methods developed in this paper will 
contribute to the understanding and a better interpretation of such experiments. 

The benefit of the quantum normal form presented in this paper is not only to give 
a firm theoretical framework for a quantum version of an activated complex but it 
moreover leads to a very efficient method for computing quantum reaction rates and 
the associated resonances. In fact, the quantum normal form computation of reaction 
probabilities and resonances is highly promising since it opens the way to study high 
dimensional systems for which other techniques based on the ab initio solution of 
the quantum scattering problem like the complex dilation method |Sim79l IRei821 
IMoi98j or the utilization of an absorbing potential [NMOlj do not seem feasible. We 
mention, that in order to compute resonances from the complex dilation method that 
are sufficiently accurate to facilitate a comparison with our quantum normal form 
computations for the three-degree-of-freedom example studied in this paper we had 
to diagonalise matrices of size 2 500 x 2 500, and this way we reached the limits of 
our numerical computation capabilities. Furthermore, the complex dilation method 
requires the "tuning" of the scaling angle which is not straightforward but has to 
be worked out by repeating the numerical computation for different scaling angles. 
In contrast to this, the quantum normal form computation can be implemented in a 
similarly transparent and efficient way as the classical normal form. The quantum 
normal form then gives an explicit formula for the resonances from which they can be 
computed directly by inserting the corresponding quantum numbers. In particular, 
this leads to a direct assignment of the resonances which one cannot obtain from the 



ab initio methods mentioned above. 

We used the Weyl calculus as a tool to systematically study several further aspects 
of the quantum-classical correspondence. One such aspect is the relation between the 
quantum mechanics of reactions to the phase space structures that control classical 
reaction dynamics. We showed that the scattering wavefunctions are concentrated on 
those Lagrangian manifolds foliating the neighbourhood of a saddle whose toroidal 
base fulfill Bohr-Sommerfeld quantisation conditions. The location of such a 'quan- 
tised' Lagrangian manifold relative to the NHIM's stable and unstable manifold, i.e. 
the question of whether the classical trajectories on such a Lagrangian manifold are 
reactive or nonreactive, determines whether the scattering wavefunction corresponds 
to an open or a closed transmission channel. In fact, the cumulative reaction proba- 
bility can be interpreted as a counting function of the number of open transmission 
channels (i.e., the number of quantised Lagrangian manifolds in the reactive volume 
of phase space) at a given energy. We showed that the Weyl approximation of this 
number is obtained from dividing the phase space volume of the invariant subsys- 
tem representing the activated complex enclosed by the NHIM of the given energy 
by elementary quantum cells, i.e. quantum cells with sidelength given by Planck's 
constant. 

We moreover showed that the resonance states can be viewed to be localised on 
Lagrangian manifolds for which in addition to the Bohr-Sommerfeld quantisation of 
the toroidal base the remaining degree of freedom fulfills a complex Bohr-Sommerfeld 
quantisation condition. These complex Lagrangian manifolds project to the NHIM 
and its unstable manifolds (the direction of the decay of the resonance states) in real 
phase space. 

Most of the theory discussed in this paper, both classically and quantum mechan- 
ically, is local in nature. In fact, the flux in the classical case, and the cumulative 
reaction probability and the associated resonances in the quantum case only require 
local information derived from properties of the Hamilton function or operator, re- 
spectively, in the neighbourhood of the saddle point. This information can therefore 
be extracted from the classical and quantum normal form. Some of the classical 
phase space structures in the neighbourhood of the saddle where they are accurately 
described by the normal form are non-local in nature. This concerns the stable and 
unstable manifolds and the Lagrangian manifolds mentioned above. In fact they can 
extend to regions far away from saddle point. This 'global' information is important 
for the study of state specific reactivity and the control of reactions. Since these 
phase space structures are invariant manifolds and hence consist of trajectories they 
can be obtained from "growing" them out of the neighbourhood described by the 
(classical) normal form by integrating the equations of motion generated by the orig- 
inal Hamilton function. For the classical case we used this as we mentioned in this 
paper to develop an efficient procedure to determine, e.g., the volume of reactive 
initial conditions in a system. For the quantum case we mention the recent work 
by Creagh [Cre041 ICreOSj who developed a semiclassical theory of a reaction opera- 
tor from a kind of normal form expansion about what we defined as the dynamical 
reaction path in this paper. Our own future work will follow similar ideas by extend- 
ing the Bohr-Sommerfeld quantised Lagrangian manifolds that carry the scattering 
wavefunctions to the Lagrangian structures associated with the asymptotic states of 
reactants and products. The goal is to develop an efficient semiclassical procedure to 
compute full scattering matrices. This would not only allow one to compute state- 
specific reactivities but also give a clearer idea of how the quantum signatures of the 
activated complex are manifested in scattering experiments. 



Appendix 



A Proof of Lemma [9 



We here provide a short sketch of the proof of Lemma O 



Proof. Let A G 5^(11'^ x R'^) and A' be the symbol of e^^^^^^ Op[^]e-i OpI'^l with 
W G W^j^.i^^. We need to show that A' G ^^(R'^ x R'^). To this end define for s > 0, 

W'(E'^) := {ip G l2(E'^) : Op[B]V' G l2(]R'^) V5 G Wqi„withO < s' < s} , (A.l) 

and let W~''(R'^) denote the dual of ?^''(R'^). Then a variant of the usual Deals 
characterisation of pseudo-differential operators gives that A G SniR"^ x E'^) if and 
only if for all s,s' G Z, 

Op[A] : W'^'(R'^) ^ W^(R'^) . (A.2) 

This follows because (IA.2P implies that for any Bj G Wqm, j = 1, • • • , with Sj > 0, 
we have 

[Op[BN],[Op[BM-i]r-- ,[0p[B,],0p[A]]]] : L^R'') ^ L\R'') , (A.3) 

and this implies that A G SniM.'^ x R'^) (see [DS99] ). 

For a real valued G W^^^.i^^, we define U{e) := e-s^°P[^l Then U{e) : 
L^(R'^) L^(R'^) since U{e) is unitary. Moreover, we have that 

f/(e) : n'(R'^) n'iR'^) (A.4) 

for all s. To see this, let V S W*(R'^). Then we have to show that Op[-B][/(e)V' G 
L2(R'^) for all B G with s' < s. But Op[S]{7(e) = U{e)U{-e) Op[B]U{e) and 

[/(-6) Op[B]U{e) - Op[B] =1^-^ {u{-e') Op[B]U{e')^ de' 

= 1^ Ui-e')^[Op[W],Op[B]]U{e') de' . 

Hence 

Op[B]U{e) = U{e)(^Op[B] + iy{-e')^[Op[W],Op[B]]U{e') de') . (A.6) 

Since W is localised, the commutator ^[Op[Ty], Op[i3]] is a bounded operator, and 
therefore Op[B]U{e)'tp G L^(R'^). 

By (IXil) we see then that if Op[A] satisfies (1X21) then U{-e) Op[A]f7(e) satisfies 
(IX2|) . too, and therefore A' G Sni'R'^ x R'^). □ 



B Symbols of the quantum normal forms of 
the systems studied in Section [7] 



Table II: Nonvanishing coefficients of the symbol HQj^p{h,x,^) = 
Ea+/3+27<io h{a,p,'y)x'^i^K^ of the 10*^ Older quantum normal form of the one DoFEckart 
barrier with the potential (17.21) studied in Sec. 17.11 Recall that the nonvanishing terms in 
the normal form have a = (3. 
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Table III: Nonvanishing coefficients of the symbol -f^QNF ~ 

X]|q|+|/3|+27<io ^(a:/3.7)'^i^'^2^^f^^2^^''^ ^^^"^ ordcr quautum normal form of the 

coupled two DoFEckart-Morse system defined in Equation (17.151) in Sec. 17.21 Recall that 
the nonvanishing terms in the normal form have a = [3. 
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Table IV: Nonvanishing coefficients of the symbol -^qnf ~ 

T.\a\+\i3\+2y<ioh»,l3r()^i'xTxT^i'^2^^i'f>'^ of the coupled 3 DoFEckart-Morse-Morse 

system defined in Equation fl7.38p in Sec. 17.31 Recall that the nonvanishing terms in the 
normal form have a = (3. 
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C Computation of Quantum Resonances from 
the Complex Dilation Method 



We here provide some details on the complex dilation method jSim79l IRei82t IMoi98| 
that we used to numericahy compute the quantum resonances of the 2 DoFcoupled 
Eckart-Morse system in Sec. 17.21 and the 3 DoFcoupled Eckart-Morse-Morse system 
in Sec. 17. 3[ We illustrate the method for the 2 DoFsystem. The generalisation to 3 
DoFis straightforward. 

Let H = Op[H] be the Weyl quantisation of the Hamilton function H defined in 
(j7.15|) . For an angle a > 0, we define the scaled operator that is obtained from 
the operator H by substituting for the coordinate x the scaled coordinate exp(ia)x, 
i.e., 

= + ^)+ ^E(e'°x) + Vuiy) - e^'e---^ . (C.l) 

For a 7^ this operator is no longer Hermitian. The effect of the complex scaling 
is that, for suitable a > 0, the generalised eigenfunctions of H that correspond to 
resonances become square-integrable after the substitution x exp(ia)x, i.e., they 
become genuine elements of the Hilbert space L^(Il^). The resonances are then 
given by the eigenvalues of the operator which can be computed from a standard 
variational principle using a finite matrix representation in which H°' is expanded in 
terms of some truncated basis set. 

We choose the basis set given by the product states \ndv,nM) '■= \ndv) {nu), 
where, using the Dirac notation, the states |ndv) and {nu) with quantum numbers 
njv and tt-m form ID basis states in the directions of x and y, respectively. For the 
y-direction, we choose the eigenstates that correspond to the discrete part of the 
spectrum of the ID Morse oscillator Hm '■= —{h^/2m)dy + Vyi{y). The quantum 
number riM then runs from to nMmax — !> where 

\/2mD(; 1 

"-Mmax = [ 7 h -J (C.2) 

aun 2 

is the number of bound states of the ID Morse oscillator. The matrix with elements 
{nM\HM\n'^) is then diagonal with the Morse oscillator energies on the diagonal, i.e.. 



a'Lh'^ ( 1 ^2^^ 



- , a? f 



2 



2 auh 

(C.3) 

In order to compute the matrix elements of H°' with respect to the product states 
we also need the matrix elements {nulPyliT''^) , where py is the momentum operator 
Py = —ihdy. For the elements above the diagonal, we get (see, e.g., |vJBv85j ) 



1 /2 



(C.4) 



where 



= 2/5 - 2nM - 1 , P = . (C.5) 

The diagonal elements vanish and the elements below the diagonal can be obtained 
from the elements above the diagonal, 

{nM\Py\nM) = 0, {nuiPyWu) = -(?^M|Py|?^M> • (C.6) 



In x-direction we choose a so called discrete value representation [LHL85j which 
consists of a basis set |ndv)i S for which the wave functions (x|ndv) are 
localised in space on a discrete grid. Concretely, we choose the "sine" functions 

/ I \ /^ sin(^(x -ndvAx)) 

X ndv = V Ax — , (C.7) 

7r(x — ndvAx) 

where Ax is a positive constant (the grid spacing). The states l^dv) are normalised 
and orthogonal. The matrix elements of the kinetic energy operator pi. /{2m) are 
easily worked out to give 

^^2 ( l^ll^ n , — n' 

\ \„> \ — ) 6mAa;2 ^ ) 'Mv — "dv (n fi\ 

(^dvl^lndv) - < / ^Nndv-n^,, ^ UA ^n' • ^ ' 

Similarly, we get for the above-diagonal matrix elements of the momentum operator 
Px in this representation 

= (-l)<v-"dvi ^ ^ > . (C.9) 

v^dv "-dvj^x 

The diagonal elements vanish and the elements below the diagonal can be obtained 
from the elements above the diagonal, 

(ndvlPxI^dv) = 0, (n'dvlPzl^-dv) = -(n-dvlPrrl^-dv) • (C.IO) 

The matrix elements of the potential V^, or more precisely the complexified potential 
V-^{x) = VE(exp(ia) x), have to be computed from numerical quadrature. 

Using the results above, we find that the matrix elements H? , , \ := 

(tt-m, '^dvl-ff'^l^dv) '^m) of the full operator iJ" are given by 

+ EuinM) ^ndv^^ Kmu'-^ + e"'" (ndvlPxIn'd^) (nM|Py|f^M> • 

(C.ll) 

In our numerical study of the 2 DoFsystem we chose hm £ {0, . . . , 13} (for our 
choice of parameters in Sec. 17.21 the Morse oscillator has 14 bound states), ridv £ 
{-50, . . . , 50}, and A^ = 0.1. This led to a matrix of size 1414 x 1414. In our 
numerical study of the 3 DoFsystem we chose ndv G {— 25, . . . , 25}, A^; = 0.16, 
'^M;2 G {0, . . . , 6} and nM;3 S {0, . . . , 6} (for our choice of parameters in Sec. 17.31 the 
ID Morse oscillators have 14 and 17 bound states, respecively). This led to a matrix 
of size 2499 x 2499. We computed the eigenvalues of these matrices using the function 
eigs in Matlab. For both systems we chose the scaling angle to be a = 1.2. 
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